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ABSTRACT 


HIGH-RESOLUTION LASER ABSORPTION 
SPECTROSCOPY OF OZONE NEAR 1129.4 CM-1 


Lawrence N. Majorana 
Old Lominion University, 1980 
Director: Dr. Gary Copeland 


A Beer's Law experiment was performed with a tunable 
diode laser to determine self broadened line shape 
parameters of one infrared absorption ozone line in the <0 i 
band for ten pressures from 0.26 to 6.29 Torr at 285K. SO 2 
line positions were used for wavelength calibration. Line 
shapes were iteratively fitted to the Voigt function at a 
Doppler width of 29.54 MHz (HWHM) resulting in values for 
the integrated line strength, (S), of 
(0.144 +/- 0.007 ) X 10"^° cm-l/molecule cm-2, line center 
frequency, (Vg) , of 1129.426 cm-1 anc the Lorentzian 
contributions to halfwidth, (ol) . A linear least squares fit 
of as a function of pressure yielded a zero intercept of 
15.27 +/- 0.29 MHz ( rho = 0.99) and a broadening parameter, 
( olq ), of 5.71 +/- 0.29 MHz/Torr. This results in a line 
width (FWHM;i of 0.144 +/- .007 cm-1 at 760 Torr and 285K. 
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INTRODUCTION 




I. 


PURPOSE 

With the advent o£ semi-conductor lasers and their 
recent application to spectroscopy [ 2 , 16, 17, 22] » there has 
been a revolution in the study o£ in£rared spectra. Use o£ 
the tunable diode laser (TDL) has permitted ultra high 
resolution (< 0.0001 cm- 1 ) studies to be per£ormed £or 
several atmospheric constituent molecules o£ high current 
interest (CIO, HNO 3 , SOj , CF 2 CI 2 , etc.) [2,17,37]. This 
study is concerned with the experimental determination o£ 
pressure broadening characteristics o£ a single ozone 
absorption line near 1129.4 cm-1 , utilizing a TDL as a 
radiation source. 

The shape o£ any spectral line contains in£ormation 
which is indicative o£ the conditions existing in the 
vicinity where the absorption or emission o£ radiation takes 
place. Consequently, spectral analysis o£ the £actors which 
contribute to the £inite width, strength, and position o£ an 
observed line leads to a determination o£ temperatures, 
presb’.’re, concentrations, and species o£ the source 
moleculss, even while existing in extraordinarily small 
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amounts or at extreme distances. Remote sensing techniques 
such as heterodyne radiometry and differential lidar have 
been developed [16# 22] in the application of calculating 
atmospheric transmittance. These systems are designed to 
retrive high resolution spectral line profiles# from which 
by means of analytical inversion methods# it is possible to 
infer concentration and temperature profiles of atmospheric 
constituents. Most of these techniques operate in spectral 
regions where the molecular species of interest exhibits a 
strong absorption coefficient and absorption by interfering 
species is minimal. In order to produce accurate results# a 
precise knowledge of spectral line parameters is required. 

A high resolution ground observation infrared Fourier 
Transform spectrogram of solar spectra around 8.85 microns 
by Goldman et.al. [12] is shown in figure 1. An indicator 
marks the ozone vibration rotational transition 
(vl#v2,v3#J#KA#KC) (1#0#0#31#1#31 «- 0 , 0 , 0 # 32 # 0 # 32 ) centered 
at 1129.426 cm-1 (8.85406 microns). This line is a good 
candidate for use in the remote sensing of ozone since it 
meets the criteria with respect to windows# strength and 
favorable ground state. 
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The subject of this thesis is an experimental 
investigation of this line which results in a determination 
of the absorption coefficients from which "independent” 
values of line strength, Lorentz self pressure broadened 
halfwidth, Doppler width, and ozone > ozone collision cross 
section are found. In addition we will explore some of the 
physical phenomena that contribute to spectral line 
broadening and the experimental method used to measure 
spectral line parameters. 


REVIEW OF RESEARCH 

A survey of relevent ozone line parameters available to 
date (mid 1980} is presented in Table 1. [15]. Column 1 of 
the table lists the author; column 2 the spectral range; 
column 3 whether band model parameter or single line 
parameter; columns 5,6,7, and 8 list values for half-width 
at half absorption maximum, a, for ozone-ozone, ozone-oxygen, 
ozone-nitrogen and ozone^^air broadening, respectively. 
Column 9 indicates whether it was an experimental (e) or 
theoretical (t) determination. 
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From an ex^unlaation of table 1, it is apparent that: 

1. ) Only 2 experimental determinations of O 3 -O 3 broadening, 

have been reported, both at microwave wavelenghts. 

2. ) No O 3 -N 2 measurements have been reported 

in any spectral region. 

3. ) No direct measurements in the infrared have been made. 

The early low resolution study by Walshaw 1954 [45] 

estimated average air broadened linewidths in the v3 
fundamental band to be .078 cm-latm-1 (HWHM) at 20°C. This 
"independent” value has been used as a standard, implicitly 
or explicitly, by many later workers. 

The millimeter wave spectrum of ozone has since been 
observed by Lichtenstein et. al. 1971 [29] and produced 
the only experimental measurements of individual linewi hs 
until this year (1980). They reported self-broadened 
linewidths for two transitions at 118.4 GHz (J,Ka,Kc) « 
(1,1,1 0,0,0) and 110.8 GHz (6,1,5 4- 6,0,6) yielding an 

average halfwidth of 4.63 +/- .24 MHz/Torr (.J17 cm-l/atm-1) 
(HWHM). The exact experimental situation is unclear from 
their paper, however, and thus leaves some question as to 
the purity and pressure of the ozone in their work. 

Tejwani and Yeung 1975 [41] present a theoretical 

calculation of O 3 -O 3 , O 3 -O 2 , O 3 -N 2 and 03 -air linewidths 

based on Tsao-Curnutte-Anderson theory for both A and B type 
bands. The average self broadened HWHM value for type A and 
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B bands was found to be .1128 cm-1 atm-1 and .1104 cm-1 
atni-1 respectively. From their tabulation it is seen that " 
a O 3 -O 3 show quite appreciable variation with respect to 
quantum numbers.” The broadening parameter value obtained by 
this calculation depends in part on the quadrupole moment 
adopted for O 3 which is not well known. Their method was to 
use a value which gave the closest fit to the experimental 
data of Lichtenstein et. al. [29]. Not surprisingly 
therefore, the data agrees well with Lichtenstein's and so 
the result is not independent. It is the Tejwani- Yeung 
average halfwidth value that is currently listed on the 
AFCRL line parameters compilation [32]. 

Recently, Monnanteuil 1980 [35] has reported privately 
unpublished values for O3O3 of 5 microwave transitions 
listed in table 1. Details of the study, however, are not 
available at this time. 


II. THEORETICAL FORMULATION 


INFRARED SPECTROSCOPY - DEFIKl i JS 

The theory and application of direct absorption 
spe-'troscopy is discussed in many references [6,14,18,42,43] 
to which the reader is referred for detailed development. 
However, in order to introduce and define some important 
concepts used here, the spectroscopic procedure and theory 
are briefly outlined below. 


1. Transmittance 

Basically, as radiation of intensity Io(v), incident on 
a transparent cell containing a gaseous sample of the 
molecule of interest is "tuned" in frequency, absorption of 
the radiation will occur at discrete frequencies which give 
rise to the molecular spectrum. At these frequencies the 
portion of the radiation absorbed is governed by the 
strength of the particular line. The ratio of incident 
intensity, Io(v), to the intensity that passed through the 
cell, I(\»), is defined as the transmittance (T). A typical 
plot of the transmission function over a single line may 
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appear as in our result (figure 11). The line center (mo) 
is the point in frequency where the transmittance is least 
(maximum absorption). 


2. Absorption Coeffici<«nt 

The transmittance is proportional to the concentration 
of molecules in the cell and is given by Beers' law 
(sometimes r^fered to as the Lambert - Bouguer law) given 
by: 


T(v) ■ I(v)/Io(v) » exp(-k(v)lp) (1) 

where 1 is the path length of absorption, p is the partial 
pressure, of the absorbing gas. and )c(v) is the absorption 
coefficient evaluated at the frequency, a . Figure 12 is a 
plot of the absorption coefficient which gives rise to the 
transmittance shown in figure 11. Commonly, a parameter 
used to quantify this profile is the full width at half 
maximum intensity (FWHM) designated by a . This parameter 
has special significance with respect to the lifetime of the 
excited state. Another parameter used to describe the 
profile is the integrated line strength and is defined as 
the integral of the absorption coefficient with respect to 
frequency. 



As a matter of reference, frequency or wavelength is 
customarily given in units that are contingent on the 


spectral region of interest. In the infrared, 

spectroscopists often use the wavenumber ( v' ) , rather than 

wavelengths. Wavenumber is the number of waves per unit 

length, i.e. the reciprocal of wavelength. The units are 

expressed in reciprocal centimeters (cm-1), while 

-20 

wavelengths are in units of microns (u) where lu « 1 xlO 
cm • 10,000 cm-1. Another name for the cm-1 is the Rayser, 
the relation being, 1 millikayser (rak) ■ 1 xlO"^ cro-1 * 30 
MHZ. 

THE OZONE MOLECULE 

Ozone is a minor constituent of the earth's atmosphere 
existing in concentrations on the order of a few ppm in a 
tenuous layer 5 to 10 km thick at an altitude of about 30 Km 
in the stratosphere. It plays an extremely important role 
in radiative transfer and in the photo-chemical cycles of 
the atmosphere. The molecule has several absorption bands 
in the ultraviolet (xOOOnm) that are responsible for 
filtering approximately 99% of the radiation that would 
otherwise reach the earth's surface. Most of this is due to 
the Hartley absorption band which consists of a strong 




continuum between 300 and 220 nm with a very high and almost 
symmetrical peak near 255 nm [19]. 

This triatomic, non-linear molecule is a slight 
asymmetric top (lai^Ibffclc) with C 2 v symmetry and an apex 
angle of 116.8°[44). Its symmetry and three fundamental 
vibrational modes are illustrated in figure 2. 




Figure 2. Normal vibrational modes of the ozone molecule. 

The resulting spectra of ozone as any asymmetric 
molecule iu fairly complex. Ozone is a slight to weak 
dipole having a dipole moment of 0.53 debye (0.53 xlO"^ ® 
esu) [44]. From microwave studies [36] fundamental 
rotational constants are determined to be; 

A « 106,534.74 MHz 
B = 13,348.95 MHZ 
C = 11,834.30 MHZ 

The rotational constants are defined as A^h/Sir^Ig, 
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Where h ie Planck's constant and 1 ^ , 
Ig, are the inertial moments about the corresponding 
axes. Prom these constants, vibrational rotational band 
centers have been designated by Kaplan et.al. [28] as: 

vi ■ 1110 cm-1 
V 2 ■ 70S cm-1 

V 3 ■ 1042.2 cm-1 

These moments of inertia yield a value of -0.97 for Ray's 
asymmetry pareuneter (K), defined as: 

(2B - A - C) 

^ T‘=“5 

Thus, for the ground state, this molecule is very nearly a 
prolate symmetric top. 

The small binding energy (Dissociation energy of 1.04 
eV [19]) accounts for the extreme reactivity so 
characteristic of this most oxidizing of all molecules. 
This fact has hindered extensive study of high pressure 
ozone due to the imposing problems of purification and 
explosive hazard. 
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THEORY OF LINEWIDTHS AND SHAPES 

Th« spectral profile of figure 11 is fundamental to our 
study. It is actually a convolution of several separate 
broadening processes. They are the natural, Doppler, and 
collisional lineshapes that are observed together and are 
referred to as the Voigt shape. Each is briefly considered 
below. 


1. Natural Lineshape 

The natural or radiative lineshape is a direct 
consequence of the finite lifetime of an atomic or molecular 
energy state. It can be viewed either in a classical sense 
or quantum mechanically, both yielding the same result. 

Classically, we can picture the system as a bound 
harmonic oscillator damped by virtue of its own radiation 
[6,27]. A good approximation of the equation of motion for 
the oscillating electron is given by: 

X + YX + cOq^X » 0 (2) 


where: 


IS 


is the classical decay rate and Wq* resonant angular 
frequency. 

For Y<< (Oq , Equation (2) has solution: 


X ■ Xoexp(-Yt/2)exp(-iu»ot) (3) 

which has the form of a damped oscillation. The electric 
field can be determined through Maxwell's equations to be: 

E(t) - E(0) exp {-i(wo “ ^)t} t i 0 

^ (4) 

E(t) “0 t i 0 

which are similarly time dependent. Since the wave train is 
not infinite, it is not truly monochromatic and can 
therefore be resolved into its component frequencies by 
taking the Fourier transform: 


E(«) 



exp (iut) 3t 


(5) 


The intensity or power distribution, Ku), is proportional 
to E(o))^and so is found to be: 


I(w) 


lo 


y/2w 

(u-(Do ) ^+Y^/4 


( 6 ) 
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or in terms of frequency: 

KW - I, ^ (7) 

(v (Y /4 

where lo is the maximum intensity at the line center 
frequency wq# or vq $ and 


■ drassleal^llfatiM ‘ ®’ 

This important functional form, (6) or (7), is called a 
Lorentzian distribution and is the natural lineshape. It i.: 
illustrated in figure 3. The linewidth FWHM is given by, 

o ■ Y ■ 1/T classical • 2 itAv (9) 

Quantum mechanically the same result can be obtained in 
a simple manner by an application of the uncertainty 
principle: 


A Eat >. ft (1.:) 

at is the time uncertainty that an atom will spend in 
an excited state and is measured by the mean lifetime of 
the state. aE is the corresponding energy uncertainty of 
the level, which is. 
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E • h V 
or 

■ ‘'’ 2 * 

So tht fr*qu«ncy spread of a particular stats, say j, can 
thus bs written: 


2vTj 

which compares to equation (9). Actually one should use an 
uncertainty of v ■ v + v , but is negligible for the 
ground or metastable state since • . 

The lifetime of an excited state is related to the 
transition probabilities by. 


T ■ 1/ r Aj i 


( 12 ) 


if more than one transition is possible. It can be shown 
[42] (in MKS units) that Aj^ is proportional to by: 



3 



coc” 


vi 

h 



(13) 


where P.w ^ is the dipole moment matrix element, 
ox 

Consequently natural line width decreases rapidly in the 
infrared. A comparison of the magnitude of natural widths 
in the infrared to optical frequencies can be gained using 
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(11), (12), and (13), i.«. 


Av « v3 (14) 

Halfwidth of a transition at tan microns comparad to ona at 
say .5 microns is approx imataly four ordars of magnituda 
smallar. Natural width in the infra;:ad will normally ba 
savaral ordars of magnituda lass than Oopplar or collision 
broadanad width. 


2. Oopplar Broadening 

In the regime of low gas pressure the observed 
linewidth is due primarily to the Doppler shifting of 
emitted or absorbed frequencies due to thermal motion of the 
molecules. Molecules moving with a component of velocity 
away from the observer give rise to a red shift and 
conversly those moving towards, shift blue. Since the shift 
varies from one molecule to another the effect is not the 
same for all molecules and therefore classified as an 
inhomogeneous broadening me. an ism. 

The observed shifted angular frequency is given by, 

Oj' ■ Ug (l-^*?/c) 


( 15 ) 
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wh«r« 9 is a unit vsctor in a dirsction from ths obssrvsr to 
th« molsculs and is the emitted or absorbed frequency of 
the stationary atom. 

Statistically tne probability of an atom having a 
velocity between v and v^dv in a one dimensional coordinate 
system is given by the Maxwellian distribution^ 


P(v)dv ■ 



(16) 


where » 

M « WiAas of the atom 
T ■ absolute temperature 
)c ■ Boltzman's constant 

Let the line width parameter ^ ^ / be defined by 



(17) 


then from (15), (16), and (17) the probability of absorbi»d or 
emitted radiation having an angular frequency between <*»o' 
and Wo' + <3 wo' is given by. 


P(uo)dwo » exp [ -4 (wq - Wg)^/4^ j du' (18) 


This form of the profile is a Gaussian distribution (due to 
the Maxwellian velocity distribution), about the central 
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frequency uo: whose width (FWHM) is given by# 
aj - Min J)** - 2 ‘ 

Compared to the Lorentzian function# the Gaussian drops 
off rapidly on either side of its center# (see figure 5). 
In fact# 99.994% of the area under the curve is contained 
within the bounds of V** four linewidths. The Doppler width 
in terms of wavelength is given by 

ttjj - 7.16 X 10”^iV T/A (20) 

It is apparent that oiq varies linearly with wavelength. For 
the ozone molecule# (A « 48)# at a temperature T « 300K# one 
finds in the ultraviolet (lOnm)# the Doppler width is 
approximately 1.8 xlO'^ nm which is on the order of the 
natural width. In the infrared (-10 u) the Doppler width of 
ozone is approximately 60 MHZ# two orders of magnitude 
greater than the natural width. Observed spectra in this 
experiment is Doppler limited. 


3. Collision Broadening 


In all real physical situations of interest a molecule 
is not alone in the universe » and so will consequently be 
subjected to interactions by neighboring atoms, ions, or 
electrons, which in turn perturb the molecular energy 
levels. For a large collection of such events averaged over 
many molecules, broadening of its spectral lines is 

observed. 

The first successful qualitative collision theory was 
formulated by Lorentz [30] in seeking explanations of 
absorption experiments by Angstrom and Hallo. He envisioned 
the excited electron oscillating with a subsequent field 
f(t) and assumed that the perturbation due to the 

interaction of a colliding molecule quenched the radiating 
wave train abruptly. Disregarding the relatively small 

radiation damping, a Fourier analysis of f(t) similar to 
that for the natural broadened line yields the power 
distribution illustrated in figure 4b. For many collisions, 
averaging over all possible values of t, produces the 
Lorentzian profile: 

I(u) - Iq U/lli (21) 

(a» - + (1 /t)^ 

shown in figure 4c, where ^ is the mean time between 
collisions and the FWHM 2 / t . 
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Figure 4. Effect of abrupt cut-off of radiation, (a) shows 
the rectangular cut-off function of duration t^. (b) 
shows the frequency spread resulting from the cut-off 
of width 2it/tu. (c) shows the Lorentzian distribution 
obtained by averaging (b) over all possible values of 
tu* ^ is the mean value of tu- 
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This method of calculation is termed the impact 
approximation and assumes that the duration of r collision 
is small compared to the time between collisions. This 
result is used in many cases to describe line broadening 
below 100 Torr. Collision theory has since been extended to 
include other effects: phase changes on impact [46,47]; 
density of perturbers; long and short range forces, and 
other mechanisms of interaction. 

Compared to the Etoppler profile, the Lorentzian 
distributes more energy away from the spectral center. 
Consequently, pressure effects are dominant in the wings of 
the line. 

The mean time between collisions, T^. , in terms of the 
mean free path is found from kinetic gas theory to be. 


T^ = X/V 


( 22 ) 


The mean relative velocity, v , is given for a Boltzman 
distribution of like collision partners as: 


V = 



(23) 


Thus, the linewidth, a[_ , is given by 


L T, 


2Nav 


(24) 
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where N is the number density and a is the collision cross 
section. 


4. The Voigt Profile 

In almost any absorption experiment, the observed 
lineshape will not be a simple Lorentzian or Guassian, but a 
convolution of the two called the Voigt profile. It is 
derived by averaging the Lorentz function (equation 21) over 
the thermal distribution of equation 16, resulting in the 
probability distribution P(x,y), (normalized to unity over 
given by Armstrong ( 3 ) : 

P(x,y)=iJ^ K(x,y) (25) 

where K(x,y) is given by; 


K(x,y) 


1. f - ^ dt 

^0 n J y2 + (X - t^) 


(26) 


and is )cnown as the Voigt function. K(x,y) is a reduced 
absorption coefficient, where: 



(27) 


equals the dimensional constant of the reduced absorption 
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coefficient. S, is the integrated line strength given by: 


S 



(28) 


It is commonly expressed in units per absorbing molecule 
( cm-l/molecule cm-2). 


X = (In 2)** 

“D 


(wave number scale in units (29) 
of Doppler width HWHM) 




2kTln 2 
Mc^ / 


(Doppler half width HWHM) 


(30) 


where is the wave number at line center and v is the wave 
number at which is to be evaluated. 

The shape of this profile is determined by the ratio of 
Lorentz to Doppler widths, designated here as y, where 


y = 



(31) 


It is this parameter that will be evaluated from the 
experimental data. Unfortunately, the Voigt profile cannot 
be expressed in a concise analytic form, but it can be 
evaluated numerically. It has been tabulated by several 
authors (7,10,25,8] whose results have been utilized by 
others to fit the observed profiles. 


n- iiii^'r III -1 vlit tifttiM "fillfcii f I III I 


Figure 5 illustrates the contribution of the Doppler 
and Lorentz profiles to the Voigt for the case of 

a^/aQ > (In 2)^. A change in the ratio quantifies the 
dominance of one case over the other. Physically "pure" 
cases are not observed. In the limits of the ratio however, 
it should be noted that the shape reduces to the independent 
broadening formula; i.e. for the "pure" Doppler regime, 
as aL-40, K(x,o) exp(-x2); similarly for "pure" Lorentz 
as qiq-^0. 


P(x- 


x/y 


(V -Vq) 


OJL ^ -±. 


(V - 




+ a, 





TUNABLE DIODE LASERS 

The radiation source used in this experiment was a 
PbSnSe diode laser manufactured by Laser Analytics, Inc. 
graciously made available by the Laser Spectroscopy Branch, 
NASA Langley Research Center. Extensive use of tunable 
diode lasers (TDL) for this type of application has been 
made since the first publication in 1970 of Doppler limited 
spectra within the 10.6 micron band of SFe (20]. An 
excellent review of these devices is presented oy Butler 
[5] . 


The semiconductor laser is mounted on an oxygen-free 
copper cold head which is cooled by a closed cycle Helium 
refrigerator, manufactured by CTI - Cryogenics Helix, 
Waltham, Ma. A forward bias applied to the diode produces 
emission at a wavelength determined by the energy gap of the 
semiconductor. Total output power is typically about one 


milliwatt continuous wave. 


30 


There are two main aspects which make the application 
of this type of device to spectroscopy advantageous. First, 
their unprecedented narrow line width (<10 MHz) is orders of 
magnitude narrower than conventional prism instruments. The 
high resolution virtually eliminates instrumental width 
limitations [40] permitting essentially undistorted 
measurements of individual line profiles so that pressure 
broadening studies can be carried out by direct linewidth 
measurements. Secondly, they are tunable in a piecewise 
fashion through most of the infrared (three microns to 
thirty microns). Tuning can be accomplished in several 
ways: varying the diode temperature, applying a magnetic 
field, applying external pressure, or sweeping the current 
through the semiconductor. A practical arrangement used 
here, was to stabilize the diode temperature at 
approximately 37° K and vary the current linearly with time 
(at a sweep rate of 2.0xl0"‘* Amp/sec) from 1.0630 to 1.0839 
Amps ( = 20 ma change) with a corresponding frequency sweep 
of approximately .12 cm-1 (3.6 GHZ). 


EXPERIMENTAL ARRANGEMENT 

A schematic of the TDL spectroscopy apparatus is shown 
in figure 6. The laser's output radiation was collected by 
a germanium f/1.5 lens, chopped, and passed through a 




TDL SPECTROSCOPY FACILITY 


I.R. ABSORPTION CELL hgOlTe DETECTOR 



Figure 6. Optical setup for simultaneous measurement of 
reference gas spectrum, test gas spectrum and 
eatalon tuning curve. 
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Czerny-Turner double pass monochromator to isolate one of 
several longitudinal mod-'s normally emitted by these lasers. 
The output from the monochromator is divided by a beam 
splitter which directs 40% through a solid 5.08 cm germanium 
etalon and is focused on a LN 2 cooled HgCdTe detector [26]. 
The transmission signal through the etalon provides a 
calibration scale to convert the laser current scale into 
relative wavenumbers, providing the free spectral range 
(FSR) is some small fraction of the tuning range of the 
mode. In this case the FSR is .02445 cro-1. In addition, 
this signal may be used to monitor the laser tuning 
characteristics which are often complicated by effects such 
as mode hopping, multimode output, and non-linearities [16]. 
The remaining 60% is directed to another beam splitter which 
directs a portion of the beam through a reference cell on to 
another detector. This cell contains SO 2 at a pressure of 
one Torr as a reference gas, whose spectrum in this region 
is precisely known [21,39]. An accurate determination of 
wavelength can be made by measuring its separation from any 
of the reference lines. The remaining radiation passes 
through the 50 cm O 3 single pass absorption cell to a third 
detector. 

Windows of the absorption cell were made of KCl 50 mm 
in diameter and 6.5 ram thick. Ozone pressure is measured by 
a 0-100 Torr Datametrics Barocell pressure sensor and 
electronic manometer. 
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This particular arrangement permitted simultaneous 
measurement of the ozone spectra, the reference gas spectra, 
and the etalon curve with the same laser sweep. 

Figure 7 illustrates the signal processing path through 
a block diagram of the electronics. Output signals from the 
three detectors are passed through appropriate electronic 
filters, phase sensitive amplifiers, and then digitized 
using a 16 channel 12 bit ADC coupled to a PDP-8 
minicomputer which records the data on 7-track tape for 
future processing. Simultaneously the data is monitored by 
an oscilloscope and a 2 channel Hewlett Packard chart 
recorder. 


PROCEDURE 

Ozone was generated from pure dry oxygen using the 
silent electrical discharge method (13), in a sealed glass 
system. The gas was condensed in a trap by liquid nitrogen 
(77K) and concentrated by pymping off excess oxygen. The 
cell was previously conditioned with ozone to minimize the 
effect of reactions with the cell walls, thus reducing 
impurities, and explosion potential. 

Prior to producing ozone, a series of laser mode scans 
was performed to locate a strong, stable operational mode of 
the diode at the desired waveler»gth. The laser control 


Control y 3 I 



Figure 7. Electronics block dlagrsn, and signal processing path 
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system is designed to automatically tune the laser over a 
predetermined frequency range once initiated. 

The first scan was made with the laser's output blocked 
to determine an optical zero. The second scan is made 
through the evacuated absorption cell in order to record an 
unattenuated power curve (lo) by which the following 
transmission curves are normalized. Ozone was then 
introduced into the cell at a pressure of 6.3 Torr and the 
first pressure scan is performed. Ozone is pumped from the 
cell for each successive scan reducing the pressure by 
approximatly one Torr for each scan. This procedure was 
repeated for ten scans. 

Throughout a scan, the amplified detector Ou '^uts, the 
TDL current, and the SO 2 and O 3 pressures were fed to the 
A/D converter and transferred to tape by the PDP- 8 . 
However, due to inconsistencies in the data tapes, (poor A/D 
performance) they were not used in the final data analysis. 
Instead the X-Y recording was digitized using a Numonics 
digitizer and recorded on magnetic tapes which were read 
into and stored on a Dec-system 10 for later analysis. 


IV. DATA ANALYSIS 


A series of Fortran programs were developed to 
efficiently reduce the spectroscopic data collected through 
the procedure described in the previous section and are 
included in Appendices A and B. The flow chart of figure 8 
illustrates their sequence and purpose. The programs QUICK 
and VFIT perform the body of data analysis. They are 
interactive programs written for use on the DEC system- 10 
with tektronics graphics routines. Their important 
functions and organization are briefly described below. 


PROGRAM QUICK 

Program QUICK provides frequency calibration for the 
laser current array, and generates from the raw data an 
array of absorption coefficients from which the line 
parameters will later be determined. 

The raw digitized data is organized into separate files 
consisting of the independent and dependent variables for 
each type of scan as illustrated in each curve of figure 9. 
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Figure 8. 


Data processing flow diagram illustrating 
sequence of analysis. 
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Figure 9. Uncalibrated raw data as collected by the PDP/8 
miniconputer and recorded on seven track tape. The 
upper trace is the SO 2 reference spectrum . The middle 
traces represent the unattenuated power spectrum (Iq)/ 
and the 1129.426 cm-1 ozone absorption line (I), at 
2.03 Torr. The bottom trace is the etalon tuning curve 
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The independent variable is always the laser current, CUR, 
in arbitrary units. The dependent variables are the 
detector signals. The detector signal through the empty 
cell is called EMPTY. The detector signal recorded with an 
optical path through the etalon is stored in the array 
ETALON. The signal corresponding to the optical zero is 
called ZERO, and the absorption cell or the reference cell 
;iignals are stored in the array GAS (the difference 
designated in the data file header record). Program QUICK 
is designed to perform all the calibration functions with 
the reference gas array first, then loop back and read in 
the absorption gas array (ozone) in its place to complete 
the analysis and thus the reason the GAS file will contain 
either the reference or the ozone gas data. 

Immediately after reading each data file, subroutine 
CURFIT is called to "dress" the data before analysis. Here 
a cut off level is determined to eliminate portions of data 
with low S/N, as well as smoothing and alignment of the data 
arrays. Within this routine another subroutine, LINFIT, is 
called which by means of a linear least squares fit 
calculates and stores the slope and intercept of the current 
versus the location in the array. This provides a check on 
the linearity of laser current over time as well as an 
efficient way to store and manipulate the current array by 
reducing storage requirements from an entire array to only 
two numbers. Normalization is accomplished by dividing the 


r 
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GAS curve by the EMPTY curve by use of the subroutine 
NORMAL. 

Subroutine CNTFNO is called to locate line centers 
(center minimums for absorption gas data and maximums for 
ETALON used later). This is done by numerically 
differentiating the entire gas curve and then determining 
where the derivative array passes through zero from positive 
to negative sense. Greater accuracy in location of line 
centers is gained by interpolation between the adjacent 
crossing points. An array of line center locations is 
passed baclc to the main program. 

At this point a plot of the normalized reference data 
with marked centers can be produced, a sample of which is 
shown in figure 10 for SO 2 as the reference gas. This 
facilitates selection of a reference line for frequency 
calibration. The user is now asked to input the absolute 
frequency of the reference line. 

In preparation for conversion of the current array to 
frequency the etalon trace is considered by subroutine 
ETPREP. After first smoothing the etalon trace the maximum 
intensity positions are determined by CNTFND. The change of 
current between peaks is calculated and stored in an array 
from which the average and standard deviation are derived. 
A plot of the etalon transmissions with marked maxima and a 
histogram of the fringe widths vs. their frequency of 









S' 


42 

occurance may be generated from this subroutine. This 
process is important to detect subtle laser mode shifts and 
nonlinearities that can occur during operation of TDLs. 

Conversion of the current array into absolute frequency 
is initiated by first transforming current to relative 
frequency (subroutine RELFREQ). A conversion factor is 
found by the ratio of the free spectral range of the etalon 
(.02445 cm-1 for our case) to the average fringe width 
previously calculated. Both the current and line center 
arrays are converted by multiplying each data value by the 
conversion factor. Transformation of this array to absolute 
values (ABFREQ) is determined by the difference between the 
absolute frequency value of the reference line and its value 
in the relative array just calculated. 

A loop to the begining of QUICK reads in the first 
ozone absorption curve as GAS which is then similarly 
chopped, smoothed, aligned and normalized. Since the 
current array is identical the calibration sections of the 
program are bypassed. 

A normalized, calibrated plot of the ozone transmission 
function is generated, an example is shown in figure 11. An 
array of absorption coefficients is calculated using the 
Beer-Lambert law (equation one) and optionally appended to a 
data file for later processing by VFIT. Transmittances may 
also be appended to a separate data file. Data files 
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remaining for different pressures may be succesively opened 
by UUICK and processed using the same calibration data. 

PROGRAM VFIT 

Program VFIT iteratively fits the Voigt function, 
K(x,y) (equation 26), to the absorption coefficients k(v ) 
calculated in QUICK and produces estimators for the 
bac)cground signal, the center frequency of the absorption 
line, Lorentz and Doppler widths (HWHM), and the integrated 
line strength. 

It utilizes the data files generated by QUICK. The 
spectral line center is determined by the maximum functional 
value of a parabolic fit to the line peak. This center 
frequency is removed from all data and the Doppler width 
(HWHM) is calculated (equation 30) for ozone at 285K. Next 
the program estimates from the data a first guess at the 
spectral line width and numerically integrates the 
absorption coefficient to obtain an estimator of the line 
strength. The data transformation is performed that 
converts the absorption coefficient versus frequency value 
to the reduced absorption coeficient (equation 26) as a 
function of frequency in Doppler units (x). 
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This observed profile is to be fitted to the function 
K(x,y) - the Voigt profile with a fixed Doppler width. The 
function K(X/y) is generated for given x,y values using the 
code described by Armstrong 13 J. At line center (x«0)/ the 
value of the reduced absorption coefficient K{0,y) is 
observed. That value is compared with a table of values of 
K(0»y) generated by the function subprogram named Kl. A 
linear interpolation is performed in the table of values of 
K(0,y) to find the closest first guess for the value y 
(equation 31). 

The value of y, as determined above, is increased by 
25% and the reduced chi squared statistic is calculated for 
the data. This statistic is calculated for decreasing 
values of the parameter y in an iterative scheme until it is 
minimized. At this point the value of y, i.e. the ratio of 
the Lorentz to Doppler width is determined. A plotting 
routine is then called for displaying the data and the 
"best" fit (see figure 13). 

Line strength is calculated over a symmetric interval 
about line center, determined by the number of data points 
above the 10% level of the pealc. Wider lines are integrated 
over a proportionatly wider data interval. 

Minor functions of the program include data weighting 
options and bac)cground corrections. 


V. RESULTS AND CONCLUSIONS 


PRELIMINARY DATA SCREENING 

The transmittance of the O 3 absorption line taken at 
2.03 Torr and 285K is shown as figure 11 which is typical of 
the ten transmission curves after normalization and 
calibration. The 1129.418 cm-1 line of SO 2 
(vl,v 2 ,v 3 , j ,ka,kc) [ 1 , 0 , 0 , 33 , 3 , 31«-0 , 0 , 0 , 34 , 2 , 32 ] was 
selected as the standard reference line because of its 
proximity to the O3 line under study and the recent 
heterodyne confirmation [39] of its position at 1129.41835 
cm-1. Using the selected SO 2 reference data , the O3 line 
center occurs at 1129.4260 +/- .0003 cm-1. The absorption 
coefficient k(v) for all pressures were derived in QUICK and 
a typical example is shown in figure 12 for the same 
experimental conditions as in figure 11 . 


ANALYSIS 


The line 
function to 
An example of 


shape factors are derived by fitting a Voigt 
the absorption coefficient as described above, 
the best (lowest reduced chi squared) fit for 
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the same experimental conditions as in figure 11 are shown 


as figure 

13. 

The Doppler 

width 

(HWHM) at 285 K 

is 

calculated 

to 

be 29.54 MHz. 

This 

value was used as 

the 


fixed Gaussian portion in the Voigt profile from which the 
Lorentzian (HWHM) contributions (a|_) were determined as a 
function of pressure. These are tabulated in Table 2 along 
th the integrated line strengths. The singular behavior 
of the 1.49 Torr data warranted its elimination in the final 
analysis . 

A linear least squares fit of the halfwidth data as a 
function of O 3 pressure shown in figure 14 yields a zero 
pressure intercept of 15.27 +/- 0.29 MHz and a slope (a^) of 
5.71 +/- 0.29 Mhz/Torr (rho=0.99). This results in a 

calculated (FWHM) line width at 760 Torr and 296 K of 
0.1445 +/-.007 cm-1 including the Doppler contribution. The 
zero pressure intercept value cannot be interpreted as the 
natural line width, since it implies an upper state lifetime 
of 5.2 nanoseconds, i.e. too fast. Instead this Lorentzian 
contribution is interpreted as the convolution of the upper 
and lower state w'dths, the laser's width (estimated to be 
5-10 MHz) and, possibly, a small contribution from power 
broadening . 

The integrated line strengths (Table 2) are found to be 
independent of pressure and equal to 0.144 +/- 0.007 xl0"^° 
cm-l/molecule-cm-2 at 285 K. 
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TABLE 2 

HALF WIDTHS HWHM AND STRENGTHS VS PRESSURE AT 285 K 


Pressure 

(Torr) 

Half 

(mKayser) 

Widths 

(MHz) 

Strengths 

(cm-2Torr-l) Units 

0.26 

. 2942 

8. 819 

. 5457 

1. 542 

0 50 

. 3012 

9.060 

. 5345 

1.510 

1.00 

. 3485 

10.45 

. 5198 

1.469 

2.03 

.4616 

13. 84 

. 5167 

1.4 60 

2.47 

. 4384 

13.14 

.4888 

1.381 

3.01 

. 5467 

16. 39 

. 5126 

1.44 8 

3.90 

. 6167 

18.49 

.4674 

1.321 

4.98 

. 7728 

23.17 

.4906 

1.386 

6. 29 

. 8366 

25.08 

.49 78 

1.406 


Units » ( cm- 1/molecule cm-2) xl0"2i 
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The observed values of K(0,y) as a function of the 
parameter y (see equation 31) are plotted in figure 15 
together with the calculated values. Generally/ there is 
good agreement with the expected broadening function. At 
the lowest pressures the observed values are somewhat lower 
than expected/ this may be due to the spectral 
characteristics of the diode laser. Notice heie that the 
value of the ratio 01 |_/qij) is approximately one, indicating 
the range of data to be comfortably within the Voigt region. 

The Voigt function K(X/y)/ evaluated at line center/ 
(X»0)/ should be inversely proportional with pressure / see 
equation 26. Indeed this is the case and is shown in figure 
16. A linear least squares fit yields: 
K(0/y) « 0.7697 - 0.0351* P(torr) with rho » -0.977. 

Ozone - ozone collision cross section for this 
transition and these conditions is calculated utilizing 
equations (23)/ (24)/ and (25) with with our halfwidth value 
of 4.32 GHz (.144 cm-latm-1) and a number density/ N/ at 
atmospheric pressure of 2.6 xlO^ ^ /cm^ atm. The resulting 
value of the cross section/ a , is 1.66 xl(T^^ cm^ . 
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SUMMARY AND COMPARISONS 

The line parameters derived in this experiment are 
summarized in Table 3 along with the limited work of others. 

This short list of available comparisons reveals slight 
discrepancies especially in the broadening parameter, <»o » 
but generally there is reasonable overall agreement. 

Line center, V q, is in slight disagreement with the 
value of 1129.442 cm-1 listed on the AFCRL Line Parameters 
Compilation (1979). The compilation, however, has recently 
been revised to our value to appear on the 1980 version 
138] . 

CONCLUSIONS 


The independent values for self broadened halfwidth and 
line strength, summarized in table two, have been determined 
for the 1129.426 cm-1 absorption line of ozone in a direct 
absorbtion experiment. With this information physical 
constants may be determined. A workable scheme has been 
outlined for collection of spectroscopic data as well as 
software for data reduction and interpretation have been 
presented . 
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TABLE 3 

COMPARISON OF EXPERIMENTALY DETERMINED LINE PARAMETERS 


Investigator 


Vo 

( cm- 1 ) 

a 

(cm-atm)-l 

® 0 

( MHz/Tor r) 

S 

Units 

Ma jorana 


1129.426 

. 144 

5. 71 

. 144 

Lichtenstein 

129] 

- 

. 117 

4.63 

- 

Monnanteuil 

135] 

- 

.124 

- 

- 

AFCRL Tape 

132] 

1129.426 

. 110 

- 

.167 


Units = ( cm-l/molecule cm-2) xlO 


-‘llf i i r' ii f i - 111 
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The results are in reasonable agreement with previous 
investigators microwave data and theoretical calculations in 
the development of Armstrong, but there is a lack of 
existing infrared data available for comparison. 

This ultra-high resolution measurement will reduce 
uncertainty in interpretations of retrived ozone data when 
applied in remote sensing applications discussed, if used. 
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APPENDIX A 
PROGRAM QUICK 


This rroi,ram is designed to perform the main data 
processing tasks as described in detail in the section DATA 
ANALYSIS. This program is written in Fortran-10 and runs on 
a DecSystem-10 time shared computer under the TopslO 
operating system. Extensive use is made of the Tektronix 
PLOT10 graphics package, as well as, data handling routines 
selected from those of Bevington [4] and the International 
Mathematics and Statistical Library Version 8.0. 


PROGRAM QUICK 

C Ai-IALYSIS PROGRAM FOR SPECTRAL ABSORBTION DATA 

C G.E. COPELAND , L.N.MAJORANA 

C OLD DOMINION UNIVERSITY 

C NORFOLK, JA 

C 

C C.H.BAIR 

C NASA LANGLEY RESEARCH CENTER 

C HAMPTON, VA 

C 

DIMENSION CUR(501 ) ,EMPTY(501 ) , GAS (501 ) ,ETALON(501 ) 

1 , ZERO(600) ,IMAX(4) ,CSI(3,2) ,CENTER(2,501 ) ,DDS(501 ) 

2 , ETCUR(501 ) 

COMMON / A / CUR , EMPTY , GAS , ET ALON , Z ERO 
COMMON/B/C3I , ETCUR , IMAX , CENTER 
DOUBLE PRECISION FILEN 
C 

C INPUT THE NAME OF THE DATA FILE 

C 

43 TYPE 44 

44 FORMAT(/,1X, ’ INPUT THE NAME OF THE DATA FILE?’,$) 
ACCEPT 4 5, FILEN 

45 FORMAT (A10) 

C 

C 

OPEN(UNIT=21 ,ACCESS='SEQIN' ,FILE=FILEN) 

5 CALL DATAIN(N,KK,K2,K3,K4 ,K10,P3) 

C 

C ALL ARGUMENTS ARE CHANGED IN SUBROUTINE 

C 

C USE CURFIT TO CHOP AND ALIGN DATA FOR EACH SCAN. 

GOTO( 10 ,20,30,40) ,KK 
1 0 CONTINUE 
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C 

CASE 


20 

CASE 


30 

CASE 


40 
CASE 

50 

CALL 
c 

CPAUSE 

55 CALL NOR'.iAL(EMPTY,GAS,ETALON,IMAX) 

C THE ALLOWING SUBROUTINE DI FFERINTIATES THE GAS CURVE 
C AND FINDS THE ABSORBTION LINE CENTERS. 

CALL CNTFND(CUR,GAS,IMAX(2) .DDS.IPTS, I ) 
IF(K10.EQ.3) GO TO 60 


c 

K10 

- 0 

IS 

NOTHING 

c 

K10 

- 1 

IS 

S02 

c 

K10 

- 2 

IS 

FREON 

c 

K10 

- 3 

IS 

OZONE 


TYPE 433 

433 FORMAT (IX, ’WOULD YOU LIKE A PLOT OF NORMALIZED 

P-AW D.ArA?’$) 

ACCEPT 431 .TOPLT 
431 FORMAT (A5) 

IF (TOPLT. EQ. ’NO’) GO TO 434 

CALL NPLOT (DDS , GAS , ETALON , IMAX , I PTS) 

4 34 NT-0 

CALL ABFREQ (DAVE ,NT , IPTS . ABCAL . lABS) 

CALL ETPREP (ETALON, DAVE, DDS) 

437 CALL ABFREQ(DAVE, NT, IPTS, ABCAL, lABS) 

C NOW READ IN OZONE DATA 

IF(KIO.EQ.O) CLOSE(UNIT-21 ) 

IF(KIO.EQ.O) GO TO 43 
60 CALL EPLOT(GAS,IPTS,2) 

CALL TAU(P3,IPTS,K2,K3,K4) 

CALL TRANS (P3 , IPTS, K2 ,K3 ,K4) 

TYPE 69 

69 FORMAT(1X, 'IS THERE ANOTHER G- S CURVE?', $) 

ACCEPT 79,ANS 
79 FORMAT(A5) 

IF(ANS.EQ. ’YES') GO TO 43 


1-1 

WHERE WE CALL EMPTY 

CALL CURFIT (KK.N, CUR, EMPTY, IMAX, CS I, SIGMA, STORE, L) 
GOTO 5 
CONTINUE 
WHERE 1-2 GAS 

CALL CURFIT(KK,N,CUR,GAS , IMAX, CSI, SIGMA, STORE, L) 
IF(K10 .EQ. 3) GOTO 55 
GOTO 5 
CONTINUE 
1-3 ETALON 

CALL CURFIT(KK ,N , CUR , ETALON , IMAX , CSI , SIGMA , STORE , L) 
GOTO 5 
CONTINUE 

1-4 ZERO (NO CURFIT) 

IMAX(4)-N 

CONTINUE 

Z E ROUT ( EMPTY , GAS , ETALON , Z ERO , I MAX ) 

move above line over to do the plot 
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CL0SE(UNIT-21) 

STOP 

END 

C - 

SUBROUTINE DATAIN(N ,KK,K2 ,K' ',K10,P3) 

DIMENSION 

CUR(501 ) ,EMPTY(501 ) ,GAS(501 ) ,ETALON(501 ) ,ZERO(600) 

1 ,ZCUR(600) 

COMMON/A/CUR, EMPTY, GAS, ETALON, ZERO 
C READ AND TYPE HEADER DATA 
C 

10 READ(21 ,21 1 )KK,N,K2 ,K3 ,K4 ,K1 0 ,TAPENO ,P3 

211 FORMAT(1X,I2,1X,I4,1X,3("2,1X),I1 „1X,A5,1X,F) 

CIO 

READ(21 ,*)KK,N,K2,K3,K4 ,K5,K6,K7,Kw,K9,K10.K1 1 ,TAPENO 
C 1 ,PRES,TOTAV,P3 

IF(K10.EQ.3) GO TO 599 
IF(KK.GT.1 )GO TO 910 
599 TYPE 900 

900 FORMAT(/ ,2X, 'READING IN DATA: ' , // , 1 3X , ' KK' , 5X , 

1 • OF DATA PTS’ ,7X, -DATE ■ ,/) 

910 TYPE 100,KK,N,K2,K3,K4 

IF( N .GT. 501) TYPE 998 

IF( N .GT. 501) STOP 

998 FORMAT (• SORRY CHARLIE TOO MANY DATA 

POINTS I ') 

C THE HEADER DATA IS AS FOLLOWS; 

C 

C KK - TYPE OF SCAN; KK-1 IS AN EMPTY CELL 

C KK-2 GAS IN CELL 

C KK-3 ETALON SCAN 

C KK-4 OPTICAL ZERO 

C KK-5 USED TO INDICATE END OF PDP 

TAPE 

C N - TOTAL NUMBER OF POINTS 
C K2 - MONTH 

C K3 - DAY 

C K4 - YEAR 

C K5 - RUN 

C K6 - ESTIMATED FREON PRESSURE IN MICRONS ... EXCEPT ... FOR 
C PDP TAPES 69667 ,RTT20 ,1 21 19 ,1 1800 , AND 222 FOR WHICH 

C IT IS GIVEN IN TORR. 

C K7 * ESTIMATED FREQUENCY (IN MORE OR LESS ARBITRARY 
UNITS) 

C K8 - STARTING CURRENT 

C K9 - DIODE TEMPERATURE 

C KIO = 1 FOR S02 (REF GAS), 2 FOR FREON, 3 FOR OZONE 

C K1 1 » ESTIMATED CELL TEMPERATURE FOR ALL PDP TAPES BEFORE 
69667 

C » THERMOCOUPLE READING FOR ALL PDP TAPES AFTER 69667 


.1 T~ I II I ' VVliiifi .itiMli n-tr^— - 
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(INC) 

C TAPENO - PDP-8 TAPE NUMBER 
C PRES - ESTIMATED TOTAL PRESSURE 
C TOTAV - TEMPERATURE AVERAGED OVER 5 CHANNELS 
C P3 - TOTAL PRESSURE (TORR) 

C 

TYPE *,N 

GO TO (20, 30, 40, 50, 60) ,KK 

C20 READ(21,*)( (CUR(J) ,EMPTY(J) , J-I ,1+4) , 1-1 ,N ,5) 

C 

C 

20 CALL READDT (CUR, EMPTY, N) 

C DO 25 J-1 ,N 

C25 TYPE 90, CUR(J) , EMPTY (J) 

GO TO 60 

C30 READ(21 ,*) ((CUR(J) ,GAS(J) ,J-I,I44) ,1-1 ,N,5) 

C 

C 

30 CALL READDT(CUR,GAS,N) 

C 

GO TO 60 

C40 READ(21,*)( (CUR( J) , ETALON ( J) , J-1 , 1+4 ) , 1-1 , N , 5 ) 

C 

C 

40 CALL READDT(CUR,ETALON,N) 

C 

GO TO 60 

C READ NO MORE THA>i 600 ZERO DATA POINTS. 

50 IF(KK.EQ.4.AMD.N.GT.600) N-600 

C51 READ(21,*)( (ZCUR(J) ,ZERO(J) ,J-l,I+4) ,1-1 ,N,5) 

C 

C 

51 CALL READDT(ZCUR,ZERO,N) 

60 CONTINUE 

90 FORMAT (F,F) 

100 FORMAT(10X,I5,8X,I5,7X,3I4) 

r. 0 FORMAT(1X, ‘THE DATA HAS BEEN READ IN’) 

C999 FORr^T(10(1X,F7.4)) 

999 FORMAT(10(F/ .4 , IX)) 

C 

CHANGED ON JULY 11 FORM 999 FORMAT(10F8 .0) TO READ THE 
C DIGITAL TABLET RESULTS 

C 
C 

RETURN 

END 

C - 

SUBROUTINE CURFIT (KK, N , Y , X , IM , C , SIGMA , STORE . L) 
DIMENSION B(2) ,IMAX(3) ,SIGMA(3) ,X(1 ) ,Y(1 ) ,C(3,2) 
CHOPS AND ALIGNS THE DATA 




66 


IF(KK .NE. 1) GO TO MO 
C SMOOTH EMPTY CELL DATA ONLY 
50 CALL SMOOTH (X,N) 

XMAX-X(1 ) 

C FIND THE ACTUAL MAXIMUM DETECTOR SIGNAL, X(KMEAN) 

DO 60 1-2, N 

IF(X(I) .LT. XMAX) GO TO 60 
XMAX-X(I) ; KMEAN-I 
60 CONTINUE 

CUTOFF IS THE DETECTOR LEVEL BELOW WHICH THE LOWER 107, 

C IS DI REGARDED TO ELIMINATE THE PORTIONS OF DATA 
C WITH LOW S/N RATIO. 

CUTOFF-0 . 1 *XMAX 

C FIND THE FIRST POINT ABOVE CUTOFF (KMIN) . 

KMIN-1 
DO 70 1-1 ,N 

IF(X(I) .LT. CUTOFF) GO TO 70 
KMIN-I ; GO TO 80 
70 CONTINUE 

80 CONTINUE 

C FIND THE LAST POINT ABOVE CUTOFF (KMAX) . 

KMAX«-N 
DO 90 1-1 ,N 
K-N+1 -I 

IF(X(X) .LT. CUTOFF) GO TO 90 
KMAX-K ; GO TO 100 
90 CONTINUE 

100 CONTINUE 

CHECK THAT KMEAN IS NOT AT THE END OF THE CURRENT SCAN 
C (WITHIN 2 POINTS OF END). IF IT IS, SET IT IN THE 
MIDDLE. 

IF(KMEAN .GT. (N-2)) KMEAN- (KMAX+KMIN ) / 2 
C SMOOTH CURVES FOR ETALON AND GAS DATA 
110 CONTINUE 

IF(KK .EQ. 3) CALL SMOOTH (X,N) 

IF(KK .EQ. 2) CALL SMOOTH (X,N) 

C DEFINE KT2, THE NUMBER OF POINTS BETWEEN KMIN and KMAX + 
1 . 

KT2-KMAX-KMIN+1 

XKMIN-KMIN 

CALL LINFIT(XKMIN,Y(KMIN) ,KT2,B,SIGMA(KK) ,2) 

C STORE SLOPE AND INTERCEPT OF CURRENT VS LOCATION IN 
C ARRAY CALCULATED IN LINFIT. 

DO 190 1-1 ,2 
C(KK,1) - INTERCEPT 
C(KK,2) -SLOPE 
190 C(KK,I)=B(I) 

COUNT (KOUNT) THE NUMBER OF POINTS GREATER THAN 3SIGMA AWAY 
FROM 

C THE CALCULATED STRAIGHT LINE. (HOPEFULLY THERE ARE 
NONE) 
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C PRINT THIS VALUE AT END OF ROUTINE. 

SIG3-3 .0*SIGMA(KK) 

KOUNT-0 
DO 200 1-1 ,N 
AMP-C(KK,1 )+C(KK,2)*I 
DIF-ABS(AMP-Yd)) 

IF(DIF .GT. SIG3) KOUNT-KOUNT+1 
200 CONTINUE 

CALCULATE WHAT THE CURRENT IS FROM THE STRAIGHT LINE FIT 
C (EMPTY CELL ONLY) 

IF(KK .GT. 1) GO TO 210 
C Y - MX + B 

C CURRENT - SLOPE X ORDINATE + INTERCEPT 
YMIN - C(KK,2)*KMIN + C(KK,1) 

YMEAN - C(KK,2)*KMEAN + C(KK,1) 

YMAX - C(KK,2)*KMAX C(KK.1) 

210 CONTINUE 

C RECALCULATE KMIN and KMAX FOR ALL SCANS FROM 
C SLOPE AND INTERCEPT STRAIGHT LINE FIT FOUND IN LINFIT. 
C X-Y-B/M - KMIN 

KMIN-(YMIN-C(KK,1 ))/C(KK,2) 

KMAX-(YMAX-C(KK,1 ))/C(KK,2) 

CHECK THAT KMIN IS NOT ZERO 

IFCKMIN .LE. 0) KMIN-1 

CHECK KMAX IS AT LEAST 2 POINTS FROM THE END DATA POINT IF 
IT 

C HAPPENS TO APPEAR THERE 

IF(KMAX .GT. (N-2))KMAX»N-2 
KMEAN-(YMEAN-C(KK,1 ))/C(KK,2) 

C SHIFT AND ALIGN THE ARRAYS! 

K-0 

DO 220 I-KMIN.KMAX 
K-K+1 

IF(I .EQ.KMEAN') STORE-K+.5 
220 X(K)-X(I) 

IMAX(KK)-K 

C IMAX IS THE ARRAY POSITION OF THE MAX CURRENT USED. 

C TYPE 230, KK.N.K, KMIN. KMAX, KOUNT, STORE, 

C 1 C(KK,1 ) ,C(KK,2) ,SIGMA(KK) 

RETURN 

230 FORMAT(6I5 ,F7 .1 , 1 PEI 2 .3 ,2E1 2 .3) 

END 

C 

C SUBROUTINE SMOOTH 

C 

C PURPOSE; 

C SMOOTH A SET OF DATA POINTS BY AVERAGING ADJACENT 

C CHANNELS 

C 

C USAGE: 

C CALL SMOOTH (Y.NPTS) 
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DESCRIPTION OF PARAMETERS: 

Y-ARRAY OF DATA POINTS 
NPTS-NUMBER OF DATA POINTS 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED: 
NONE 

SUBROUTINE SMOOTH (Y.NPTS) 

DIMENSION Y(1) 

11 IMAX-NPTS-1 

YI-YO ) 

21 DO 24 1-1 ,IMAX 

YNEW-(YI+2.*Y(I)+Y(I+1 ))/4 . 

YI-Y(I) 

24 Y(I)-YNEW 

25 Y(NPTS)-(YI+3.*Y(NPTS))/4 . 

RETURN 

END 

C 

SUBROUTINE LINFIT(XI , Y ,N , B , SIGMA .NTIMES ) 
DIMENSION Y(t ) 

C THIS ROUTINE CHECKS THAT THE CURRENT IS RECORDED 

C AT A CONSTANT RATE; THAT IS, CHECKS CURRENT VS TIME 

C OR POSITION IN ARRAY, BY A LEAST SQUARES FIT. 

C GIVES THE SLOPE AND INTERCEPT AS OUTPUT RESULTS (B) . 

SIGMA-1 .OE+08 
10 M-0 

SX-0 .0 : SY-0 .0 ; SXX-0 .0 ; SXY-0 .0 
DO 30 1-1 ,N 

IFCSIGMA .GT. 1 .OE+07) GO TO 20 
F-B(1 )+B(2)*(XI-1+I) 

DIF-ABS(F-Yd)) 

IF(DIF .GT. SIG3) GO TO 30 
20 X-XI-I+I 

SX-SX+X 
SY-SY+Y(I) 

SXX-SXX+X*X 

SXY-SXY+X*Y(I) 

M-M+1 

30 CONTINUE 

B (2 ) - (M*SXY-SX*SY) / (M*SXX-SX*SX) 

B(1 )-(SY-B(2)*SX)/M 
SUM-0.0 

IFCSIGMA .GT. 1 .OE+07) GO TO 50 

DO 40 1-1 ,M 

F«B1+B2*(XI-1+I) 

DIF-ABS(F-YCI)) 

IFCDIF .GT. SIG3) GO TO 40 
F-B(1 )+B(2)*(XI-1+I) 

SUM-3UM+(F-Y(I))*(F-Y(D) 
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40 CONTINUE 

SIGMA-SQRT (SUM/M) 

ICOUNT-N-M 
GO TO 70 

50 DO 60 I-l ,N 

F-B(1 )+B(2)*(XI-UI) 

60 SUM-SUM-f(F-Y(I))*(F-Y(I)) 

SIGMA-SQRT (SUM/N) 

SIG3-3.0*SIGMA 
B1-B(l) ; B2-B(2) 

IFCSIGMA .GT. 1.0E+07) SIGMA- 1 .OE+06 
IF(NTIMES .EQ. 1) GO TO 70 
GO TO 10 
70 RETURN 

END 

C 

SUBROUTINE NORMAL (EMPTY ,GA', , ETALON , IMAX) 

C PURPOSE: 

C TO NORMALIZE GAS AND ETALON DATA 
C FOR PLOTTING NORMALIZED TRANSMITTANCE PLOTS. 

DIMENSION EMPTY ( 1 ) , GAS ( 1 ) , ETALON ( 1 ) , IM/.X ( 1 ) 
IF(IMAX(2) .NE. IMAX(3) ) TYPE 300 
300 FORMATC IN NORMAL IMAX(2) .NE. IMAX(3)') 

C 

CHECK TO SEE IMAX(1 )-IMAX(2)-IMAX(3) 

C IF NOT FORCE THEM TO BE EQUAL TO THE SMALLEST OF 

THEM 

IM-MIN0(IMAX(1 ) ,IMAX(2) ,IMAX(3) ) 

C 

C FORCE ALL TO BE IM 

IMAX(1 )-IM 
IMAX(2)-IM 
IMAX(3)-IM 
C 
C 

c 

c 

c 

DO 10 1-1 ,IMAX(2) 

10 GAS(I) - GAS (I) /EMPTY (I) 

DO 20 I-1,IMAX(3) 

ETALON (I) - ETALON (I) /EMPTY (I) 

20 CONTINUE 

C PAUSE 

RETURN 
END 

C 

SUBROUTINE CNTFND (CUR , GAS , NP'^S , DDS , IPTS , KT) 
DIMENSION 

DS(501 ) .CENTER (2. 501 ) ,BPAR(4) , C(501 ,3) , 
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1DDS(501 ) ,GAS(1 ) ,CUR(1 ) ,ETCUR(501 ) , CSI (3 , 2 ) . IMAX(4 ) 
COMMON/B/CSI , ETCUR , IMAX . CENTER 
C PURPOSE; TO FIND LINE CENTERS OF GAS DATA 
C PASSED TO ROUTINE: 

C DS - ARRAY OF POINTS 1ST DERIVATIVE SPECTRUM 

C NPTS - LENGHT OF THE ARRAY 

C KT - TYPE OF SCAN IDENTIFIER. IF KT-2 THIS ROUTINE 

C WILL FIND CENTER MAXTMUMS. KT-ANYTHING ELSE IT 

C FINDS MINIMUMS. 

C 

C PASSED BACK TO CALLING PROGRAM: 

C 

C CENTER - ARRAY OF LINE POSITIONS OF THE LINE 

CENTERS 

C 0,501) - GAS CENTERS DATA 

C (2,501) - ETALON PEAKS DATA 

C IPTS - NUMBER OF LINES - NUMBER OF ENTERIES IN 

CENTER 
C 

C (REQUIRED FOR SUPRESSION OF 2ND DIFF. CALCULATION) 

IC-501 

C 

CSAVE THE ORIF NGONAL ARRAY AND USE DDS AS A NEW CURRENT 

ARRAY 

C 

C TYPE FIRST 10 ENTRIES OF CUR GOING INTO FILL 

C DO 76 12-1 ,10 

C76 TYPE 78, (I2,CUR(I2)) 

CALL FILL(NPTS,CUR,DDS) 

C TYPE FIRST 10 ENTRIES OF NEW ARRAY DDS COMING FROM 

FILL 

C DO 77 11-1 ,10 

C77 TYPE 78 , (II ,DDS(I1 )) 

78 FORMAT (IX, 14 ,F) 

DO 601 IK-1 ,4 
601 BPAR(IK)-0.0 

C 

C SMOOTH THE DATA 1ST 

N SMOOT- 5 

DO 330 IK-1 ,NSMOOT 
330 CALL SMOOTH (GAS, NPTS) 

CALL ICSICU (DDS, GAS, NPTS, BPAR,C,IC,IER) 

CALL DCSEVU (DDS, GAS, NPTS, C,IC, DDS, DS, NPTS, DDS, 0, lER) 
C THIS LOOP FINDS WHERE LINE CENTERS OCCUR BY CHECKING 
C WHERE THE FIRST DERIVATIVE ( FOUND IN DCSEVU) PASSES 
C THROUGH ZERO FROM +TO 
DO 332 IK-1 ,NSMOOT 
332 CALL SMOOTH (DS, NPTS) 

IPTS-0 

DO 60 NPT-1 ,NPTS-1 
PT1 = DS(NPT) 


71 


PT2 - DS(NPT-*-1) 

IF(KT.EQ.2)G0 TO 10 
C MIN FINDING LOOP 

C DETERMINES WHERE THE DIFFERINTIATED ARRAY (DS) 

C PASSES THROUGH ZERO FROM + TO - . 

IF(PT1-PT2)40,40,60 

40 IF(PT1 .LT.0.0 .AND.PT2.GT.0.0) GOTO 50 

GO TO 60 

C MAX FINDING LOOP 

C DETERMINES WHERE THE DIFFERITIATED ARRAY (DS) 

C PASSES THROUGH ZERO FROM + TO - . 

10 CONTINUE 

IF(PT1-PT2)60,45,45 

45 IF(PT1 .GT.O .0. AND. PT2.lt. 0.0) GO TO 50 

GOTO 60 

CHECK THE 2ND DERIVATIVE TO SEE IF <0.0 THEN A MINIMUM FOR 
SURE 

C50 IF( DDS(NPT+1) .GE. 0.0) GOTO 60 

50 CONTINUE 

C 

IPTS-IPTS+1 
C INTEPOLATE 

C 

CENTER(KT,IPTS)-DDS(NPT)-(PT1 /(PT2-PT1 ) )*(DDS(NPT+1 ) -DDS (NPT)) 
C 

CENTER (KT,IPTS)-CUR(NPT)-(PT1 / (PT2-PT1 ) )*(CUR(NPT+1 ) -CUR(NPT)) 


60 CONTINUE 

TYPE 62,IPTS,CENTER(KT,IPTS) 

62 FORMATdX,' ',14.' CENTERS LAST LINE IS ',?,/) 

IF( IPTS .EQ. 0) TYPE 8883 

8883 FORMATC NO LINE CENTERS FOUND CONTINUUM') 

RETURN 

.END 

C 

SUBROUTINE NPLOT (CUR , GAS , ETALON , IMAX , IPTS) 

DIMENSION CUR(1 ) ,GAS(1 ) ,ETALON(1 ) ,IMAX(3) ,ISTRIN(15) 
1 ,IVERT(13) ,CENTER(2,1 ) .IXLAB1 (7) 

DATA 

ISTRIN/78, 79. 82, 77, 65, 76, 73. 90, 69, 68, 32, 80, 76, 79. 84/ 

DATA IVERT/84 ,82,65,78,83,77,73,84 ,84,65,78,67,69/ 
DATA IXLAB1 /67,85,82,82,69,78,84/ 

C PLOTS NORMALIZED DATA 
CALL INITT(120) 

AMIN-1000 . 

AMAX— 100000 . 

CALL BINITT 
CALL NPTS(IMAX(2)) 
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CALL MNMX( GAS, AMIN, AMAX) 

CALL DLIMY(0. .AMAX) 

CALL XFRM(2) 

CALL YFRM(2) 

CALL CHFXK(CUR.GAS) 

CALL DSPLAY (CUR, GAS) 

IF( IPTS .GE.1)CALL ARROW(IPTS , 1 ) 

CALL ARROW ONLY IF THERE ARE REALLY SEPERATE LINES 
CALL NPTS(IMAX(3)) 

CALL LINE(2) 

CALL CPLOT(CUR.ETALON) 

CALL MOVABS(650 ,730) 

CALL HLABELdS.ISTRIN) 

CALL MOVABS(400,50) 

CALL HLABEL(7,IXLAB1 ) 

CALL MOVABS(25,500) 

CALL VLABEL(13,IVERT) 

CALL BELL 
CALL TINPUT(ANS) 

RETURN 

END 

C 

SUBROUTINE ARROW(IPTS , KT) 

DIMENSION CSI(3,2) ,CENTER(2,501 ) ,IMAX(4) ,ETCUR(501 ) 
COMMON/B/CSI , ETCUR , IMAX , CENTER 
C PLOT A LINE AT THE POINTS WHERE WE HAVE FOUND A 

C SPECTRAL LINE BY LOOKIG AT THE 1ST DERIVATIVE 

C 

CALL NPTS(IPTS) 

YMAX- . 1 0 
YMIN-.05 
DO 10 1-1 ,IPTS 

CALL MOVEA(CENTER(KT,I) ..YMAX) 

CALL DRAWA (CENTER (KT,I) ,YMIN) 

1 0 CONTINUE 

RETURN 
END 
C 

C 

C 

SUBROUTINE FILL (NPTS , X , XNEW) 

C 

DIMENSION X(1 ) ,XNEW(1 ) 

C 

C PURPOSE; TAKE AN ARRAY (X) THAT HAS REPEATING 

SEQUENCES 

C OF NUMBERS AND CHANGE IT SUCH THAT IT IS 

MONOTON I CALLY 

C ASCENDING IN VALUE 

C 
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C RETURNED ARRAY IS KNEW 

C 

C SAVE THEOLl) ARRAY 

C 

DO 300 1-1 ,NPTS 
300 XNEW(I)-X(I) 

FS-XNEW(1 ) 

NSAME-0 

DO 20 I-2.NPTS 

IF( XNEW(I)-FS) 8,7,10 

7 NSAME-NSAME+1 

IF( I .EQ. NPTS) GOTO 14 
GOTO 20 

10 FS-XNEW(I) 

IF( NSAME .GT.1 .AND. I .LT. NPTS) GOTO 11 
IF (NSAME .EQ.O) GOTO 20 
C 

C JUST 2 NUMBERS ARE THE SAME 

XNEWCH-XNEW(I)-XNEW(I-1 ) 

XNEW(I-1 )-XNEW(I-1 )+XNEWCH/2 . 

NSAME-0 
GOTO 20 

1 1 CONTINUE 
NT-0 

DO 13 J-I-NSAME,I-1 
NT-NT+1 

XNEW(J)-XNEW(J)+NT*(XNEW(I)-XNEW(I-NSAME-1 ))/FLOAT(NSAME+1 ) 

13 CONTINUE 
NSAME-0 
GOTO 20 

8 WRITE(5,100) 

100 FORMATC ARRAY IS DECREASING IN 

VALUE********** ■ ) 

20 CONTINUE 

RETURN 
C 

C MAKE SURE LAST ONE IS DIFFERENT THAN NEXT TO LAST 

ONE 

14 NT-0 

DO 16 J-NPTS-NSAME+1 ,NPTS 
NT-NT+1 

16 XNEW(J)=XNEW(NPTS-NSAME)+NT*( 

XNEW(NPTS)-XNEW(1 ))/FLOAT(NPTS) 

RETURN 

END 

C - ------ 


SUBROUTINE ETPREP (ETALON , DAVE , DDS) 
DIMENSION 

DDS(1 ), CENTER (2, 501 ) ,D( 100) , FREQ (100) ,WTH(100) 
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1 ,CSI(3,2) ,ETCUR(501 ) ,ETAL0N(1 ) ,IMAX(4) 

COMMON /B/CS I , ETCUR , I MAX , CENTER 
C THIS ROUTINE CONSIDERS THE ETALON TRACE 

ALONE PREPARING 

C ’ IT FOR LATER CONVERSION TO RELATIVE AND ABSOLUTE 
C FREQUENCY. ETALON PEAKS AND THE CHANGE OF CURRENT 

BETWEEN 

C THEM IS FOUND THEN A SHORT STATISTICAL ANALYSIS IS 

C PERFORMED. A HISTOGRAM PLOT IS OUTPUT , USEFUL TO 

C DETERMINE MODE SHIFTS. 

C RECALCULATE ETALON CURRENT ARRAY FROM SLOPE and 

INTERCEPT. 

NPTS - IMAX(3) 

DO 5 NE-1 .NPTS 

5 ETCUR(NE) - ( (CSI (3 . 2 ) )*FLOAT(NE) )-»-CSI (3 , 1 ) 

C DO 73 1-1 ,10 

C TYPE 6,I,ETCl:F(I) 

73 CONTINUE 

C SHItT ARRAYS IF NECESSARY, SO CURRENT BEGINS AT ZERO. 
IF(ETCUR(1 ) .LT.O.O)GO TO 41 
GO TO 51 

41 SHIFT - ETCUR(l) 

DO 51 1-1 ,NPTS 
ETCUR(I) - ETCUR(I)-SHIFT 
51 CONTINUE 

C DO 7 1-1 ,10 

C TYPE 6,I,ETCUR(I) 

6 FORMAT ( 1 X, 'ETCUR( '.12 ,’ ) - ',F) 

7 CONTINUE 

C SMOOTH THE DATA 

N SMOOT - 10 
DO 330 IK-1 .NSMOOT 
330 CALL SMOOTH (ETALON, NPTS) 

C DO 9 1-1 ,10 

C TYPE 8, 1, ETALON (I) 

8 FORMAT(1X, 'ETALONC ,12, ') - ’,F) 

9 CONTINUE 

C FIND CENTERS OF ETALON PEAKS. 

CALL CNTFND (ETCUR .ETALON .NPTS , DDS , lEPTS ,2) 

C FIND CHANGE OF CURRENT BETWEEN PEAKS, STORE IN 

ARRAY DO . 

NUM-IEPTS-1 
DO 10 N-1 ,NUM 

D(N) - CENTER(2,N-M )-CENTER(2,N) 

10 CONTINUE 

C FIND AVERAGE CHANGE OF CURRENT BETWEEN PEAKS 

(DAVE) . 

DTOT =» 0. 

DO 20 M-1 ,NUM 
20 DTOT=DTOT+D(M) 

DAVE-DTOT/FLOAT(NUM) 
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TYPE 435 

435 FORMATOX, 'WOULD YOU LIKE A SHORT STATISTICAL 
PLOT 

1 OF THE ETALON TRACE? '$) 

ACCEPT 4 36, ST AT 

436 FORMAT (A5) 

IF (STAT.EQ. 'NO' ) GO TO 437 
C PLOT ETALON TRACE AND CENTERS FOUND 

CALL EPLOT (ETALON, lEPTS, 1 ) 

DO 14 N-1 ,NUM 
TYPE 15,N,D(N) 

15 F0RMAT(1X, 'D(',I2, ') - ',F) 

14 CONTINUE 

TYPE 25 DAVE 

25 FORMATdX, 'DAVE - ’ ,F) 

C FIND STANDARD DEVIATION, (SIG) 

SUM-0 . 

DO 30 I“1 ,NUM 
SUM-SUM+(D(I)-DAVE)**2 
30 CONTINUE 

SIG-SQKT(SUM/(FLOAT(NUM)-1 .) ) 

TYPE 40, SIG 

40 FORMAT(1X, 'STANDARD DEVIATION IS'.F) 

PAUSE 

C AT THIS POINT WIDTHS OF ETALON PEAKS WILL BE 

C COMPARED IN ORDER TO DETECT SUBTLE MODE SHIFTS 

C THAT MAY HAVE BEEN PREVIOUSLY OVERLOOKED. THE 

C WIDTHS WILL BE SORTED AND PLOTED VS THEIR 

C FREQUENCY OF OCCURANCE, (HISTOGRAM). 

C THE DISTRIBUTION SHOULD BE GAUSSIAN! 

C FIND MIN AND MAX ETALON PEAK WIDTHS , (DMIN and 

DMAX) . 

DMIN-D(1 ) 

DO 50 K»2,NUM 

50 IF(D(K) .LT.DMIN)DMIN-D(K) 

DMAX-D(1 ) 

DO 60 L-2,NUM 

60 IF(D(L) .GT.DMAX) DMAX-D(L) 

C WHAT IS THE RANGE IN VALUES OF D? 

RANGE - DMAX- DMIN 
C BUILD ARRAY OF SLOTS. 

XINC - RANGE / FLOAT(NUM-I) 

WTH(1 )-DMIN 
DO 70 N-2,NUM 

70 WTH(N) = WTH(N-1 )+XINC 

C SORT ETALON WIDTHS AND FIND THEIR FREQUENCY OF 

OCCURANCE, FREQ(I). 

DO 90 M-1 ,NUM 
1 = 1 

IF(D(M) .LE.WTH(D) GO TO 87 
I - 1 + 1 


85 








T 
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GO TO 85 

87 FPi:Q(I)-FREQ(I)+1 . 

90 CONTINUE 

C CALL PLOT ROUTINE TO PLOT HISTOGRAM. 

CALL HISPLT(WTH,FREQ,NUM) 

437 RETURN 

END 

C 

SUBROUTINE HISPLT(WTH , FREQ , NUM) 

DIMENSION FREQ(I) ,WTH(1) ,IYLAB(17) ,IXLAB(12) 

DATA 

IYLAB/70,82,69,81 ,32,79,70,32,79,67,67,85,82,65,78,67,69/ 
DATA IXLAB/69,84,65,76,79,78 ,32 ,87 ,73 ,68 ,84 ,72/ 

TYPE 10, NUM 

1 0 FORMAT ( ■ Fi'vEQ WIDTH NUM- ’ , I ) 

DO 20 1-1 ,NUM 

20 TYPE 21 , (FREQ(I) ,WTH(D) 

21 FORMAT(1X,F,1X,F) 

C PAUSE 

C PLOT THE HISTOGRAM 

CALL IN ITT (120) 

CALL BINITT 
CALL NPTS(NUM) 

CALL VBA:.ST(8,0,0) 

CALL XFRM(1) 

CALL YMFRM(I) 

CALL CHECK(WTH,FREQ) 

CALL DSP LAY (WTH, FREQ) 

CALL MOVABS(400,20) 

CALL HLABEL(12,IXIAB) 

CALL MOVABS(25 ,575) 

CALL VLABEL(17,IYLAB) 

CALL BELL 
CALL TIN?UT(ANS) 

RETURN 

END 

C - - 

SUBROUTINE EPLOT(ETALCN , lEPTS , ITYFG) 

DIMENSION ETCUR(501 ) ,ETALON(1 ) ,ISTRIN(1 5) ,CSI(3,2) 

1 ,IVERT(13) ,CENTER(2,501 ) ,IMAX(4) 

DIMENSION IVERT3(15) ,IETSTR(23) 

DIMENSION IXLAB1 (7) ,IXLAB3(13) 

COMMON/ B/CSI , ETCUR , IMAX , CENTER 
DATA 

ISTRIN/78, 79, 82, 77, 65, 76, 73, 90, 69, 68, 32 ,,80,76, 79, 84/ 

DATA IVERT/84, 82, 65, 78, 83, 7 7, 73, 84, 84, 65, 78, 67, 69/ 
DATA lVERT3/65,66,83,32,67,79,69,70,32, 

1 66,65,83,69,32,69/ 

DATA IXLAB1 /67,85,82,82,69,78,84/ 

DATA IXLAB3/65,66,83,32,70,82,69,81 ,85,69,78,67,89/ 
DATA 
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IETSTR/69,84.65./6,79.7o,32.65,78,68,32.80,69.65,75. 

1 32,67.69,78.84,69,82.93/ 

C PLOTS NORMALIZED GAS DATA 
IVALUE - IMAX(3) 

CALL INITT(120) 

CALL BINITT 
CALL NPTS (IVALUE) 

CALL XFRM(2) 

CALL YFRM(2) 

CALL CHECK (ETCUR.ETALON) 

CALL DSPLAY (ETCUR.ETALON) 

IF( lEPTS .C^.1 )CALL ARROW(IEPTS,2) 

CALL ARROW ONLY IF THERE ARE REALLY SEPERATE LINES 
CALL MOVABS(600 ,730) 

IFdTYFG.EQ.I )CALL HLABEL(23 , lETSTR) 

IFdTYFG.EQ.I .0R.ITYFG.EQ.3)G0 TO 10 
CALL HLABEL(15,ISTBIN) 

10 CALL MOVABS(25,500) 

IFdTYFG.EQ.I .OR.ITYFG.EQ.2)CALL VLABEL( 1 3 , IVERT) 
CALL MOVABS(25 .600) 

IFdTYFG .EQ. 3) CALL VLABEL ( 1 5 , IVERT3 ) 

C 

C NOW PUT IN THE X LABELS 

C 

CALL MOVABS(425,25) 

IFdTYFG.EQ.I ) CALL HLABEL( 7 , IXLABl ) 
IF(ITYFG.EQ.3.0R.ITYFG.EQ.2)CALL HLABEL( 1 3 , IXLAB3) 

C 

C 

CALL BELL 
CALL TINPUT(ANS) 

ITYFG-0 

RETURN 

END 

C 

SUBROUTINE RLFREQ(DAVE , IPTS) 

DIMENSION CENTEP.(2,501 ) ,ETCUR(501 ) , CSI (3 , 2 ) , IMAX(4 ) 
COMMON /B/CS I, ETCUR.IMAX, CENTER 
C THIS ROUTINE CONVERTS THE CURRENT ARRAY TO RELATIVE 
FREQUENCY 

NPTS - IMAX(3) 

C THEKE ARE .02445 CM-1 BETWEEN PEAKS. 

C FiND A CONVERSION FACTOR (CONF) 

CONF -.0244 5 /DAVE 

C CONVERT CURRENT TO RELATIVE FREQS. 

DO 10 J-1 ,nf: 

10 ETCUR(J) - ‘ MR(J)*CONF 

C DO 30 I-1.10C.. J 

C TYPE 20,I,ETCl’R(I) 

20 FORMATdX, ’ETCUR( ' , 12 , ' ) - ',F) 
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30 CONTINUE 

CONVERT GAS CENTER ARRAY TO REL FREQS. 

DO 40 K-1 .IPTS 

CENTERd.K) - (CENTERO ,K))*CONF 
40 CONTINUE 

RETURN 
END 

SUBROUTINE ABFREQ(DAVE ,NT , IPTS , ABCAL, lABS) 

DIMENSION CSI(3,2) ,ETCUR(501 ) ,IMAX(4) ,CENTER(2 , 5C1 ) 
COMMON / B / C S I , ET CU R , I MAX . C ENTE R 
ABFREQ CONVERTS THE CURRENT ARRAY FIRST TO RELATIVE 
FREQUENCIES (CABLING RELFREQ) ,THFN TO ABSOLUTE 

IT ASKS WHICH LINE OF THE REFERENCE SPECTRUM IS TO 

USED FOR FREQUENCY CALIBRATION AND ITS KNOWN CENTER 

IF(NT.EQ.1 )GOTO 55 
TYPE 9 

FORMAT(/ ,3X, 'CENTER 6X CURRENT POSITION’,/) 

DO 20 I-l ,IPTS 
TYPE 10, I, CENTER (1 ,1) 

F0RMAT(1X, 'CENTERC ,13, ') - ’ , F) 

CONTINUE 

C INPUT ABSOLUTE LINE FOR CALIBRATION. 

TYPE 30 

30 FORMAT(/ , IX, 'INPUT WHICH CENTER IS TO BE USED 

1 FOR ABSOLUTE FREQ. CALIBRATION? ' $) 

ACCEPT *,IABS 
TYPE 35 

35 FORMATC/ , IX, 'WHAT IS ITS ABSOLUTE FREQUENCY ? ’ $) 

ACCEPT *, ABCAL 
NT-1 
RETURN 

CHANGE CURRENT ARRAY TO RELATIVL FREQUENCY. 

55 CONTINUE 

CALL RLFREQC DAVE, IPTS) 

DIFF - ABCAL - CENTER(1 ,IABS) 

CHANGE RELATIVE ARRAY TO ABSOLUTE FREQS 
DO 60 L-1 ,IMAX(3) 

ETCUR(L)-ETCUR(L)+DIFF 
60 CONTINUE 

C DO 70 1-1 ,10 

C TYPE 75,I,ETCUR(I) 

75 FORMAT(1X, 'ETCUR( ' , 12 , ' ) = ’,F) 

70 CONTINUE 

RETURN 
END 

C 


C 

FREQS. 

C 

BE 

C 

FREQ. 
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10 

20 
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SUBROUTINE TAU(P3 , IPTS , K2 , K3 , K4 ) 

DIMENSION COF(501 ) ,CUR(501 ) ,EMPTY(501 ) ,GAS(501 ) 
DIMENSION CSK3.2) ,ETCUR(501 ) ,IMAX(4) 

COMMON /A/ CUR , EMPTY , CAS 
COMMON / B / CSI , ETCUR , IMAX 
DOUBLE PRECISION FILEA 
C PURPOSE; 

C TAU CALCULATES AN ARRAY OF ABSORBTION COEFFICIENTS 
C USING THE BEER- LAMBERT LAW: 

TRANSMITTANCE-EXP (-COF) (XL) (P3) 

C THEN APPENDS THIS DATA TO A FILE (OPTIONAL). 

C MAJOR VARIABLES: 

C XL - CELL LENGTH ^CM) 

C COF - ABSORBTION COEFFICIENT 

C EMPTY (K) /GAS (K) - TRANSMITTANCE 

C P3 - GAS PRESSURE IN CELL 

C IMAX (2) - OF DATA PTS IN GAS CURVE. (FROM CURFIT) 
CELL LENGTH IS SET HERE 
XL- 50.0 

CALCULATE ABSORBTION COEFFICIENTS. 

DO 10 K-1 ,IMAX(2) 

COF(K)-(-ALOG(GAS(K) ) / (P3*XL) ) 

10 CONTINUE 

CALL EPLOT(COF,IPTS,3) 

C APPENDING ROUTINE. 

TYPE 61 

61 FORMAT(1X, 'DO YOU WANT ABS COFS APPENDED TO A DATA 

FILE? ’ $) 

ACCEPT 62.ANS 

62 FORMAT (A5) 

IF(ANS.EQ. ’NC*) GO TO 75 
TYPE 63 

63 FORMAT(1X, ' FILE NAME? ■$) 

ACCEPT 622, FILEA 

622 FORMAT(AIO) 

OPEN(UNIT-21 .ACCESS- -APPEND’ .FILE-FILEA) 

C 

C ONLY WRITE OUT MULTIPLES OF 5 DATA PAIRS 

C 

ITES-IMAX(2)/5 

IT-5*ITES 

C TYPE *,IMAX(2) .ITES.IT 

WRITE(21 ,*)IT,K2,K3,K4,P3 
C APPEND TRANSMITENCES TO DATA FILE 
C 

WRITE (21 ,999) f (ETCURd ) , COF(I ) , I-J , J+4 ) J = 1 ,IT,5) 
999 FORMAT(5(1X,F.1X,F)) 

75 CONTINUE 

RETURN 
END 
C 
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C 


C 

SUBROUTINE READDT(X, Y ,N) 

DIMENSION X(1),Y(1) 

READ(21 , END-22, 334H (X(I) , Y(I) , I-J , J+4 ) . J-1 ,N,5) 

334 FORMAT(10(F7.4,1X) ) 

C WRITE(5,*)( (X(I),Y(I),I-J.J+4),J-1 ,N,5) 

RETURN 

22 CONTINUE 

TYPE 33 

33 FORMATC RAN OUT OF DATA IN READDT--- RECAP 

FOLLOWS ' ) 

TYPE 34 ,N 

34 FORMATC NUMBER OF DATA POINTS 

WRITE(5,*)( (X(I),Y(I),I-J.J+4) ,J-1 ,N,5) 

RETURN 

END 

C 

SUBROUTINE TRANS (P3 , IPTS , K2 , K3 , K4 ) 

DIMENSION CUR(501 ) ,EMPTY(501 ) ,GAS(501 ) ,TRA(501 ) 
DIMENSION CSK3.2) ,ETCUR(501 ) ,IMAX(4) 

COMMON / A/ CUR , EMPTY , GAS 
COMMON / B / CS I , ET CU R , I M.\X 
DOUBLE PRECISION FILEA 
C PURPOSE; 

C APPENDS TRANSMITTANCE DATA TO A FILE (OPTIONAL) . 

C MAJOR VARIABLES; 

C TRANSMITTAT^CE (TRA) -GAS (K) /EMPTY (K) (CALCULATED IN 

NORMAL) 

C IMAX(2) = OF DATA PTS IN GAS CURVE. (FROM CL'RFIT) 
TYPE *,IMAX 

C DO 10 K-1 ,IMAX(2) 

C TRA(K) = GAS (K) /EMPTY (K) 

10 CONTINUE 

C APPENDING ROUTINE. 

TYPE 61 

61 FORMAT(1X,' WANT TRANSMITTANCES APPENDED TO A DATA 

FILE? '$) 

ACCEPT 62.ANS 

62 FORMAT (A5) 

IF(ANS.EQ. 'NO') GO TO 75 
TYPE 63 

63 FORMAT(1X, 'FILE NAME? '$) 

ACCEPT 622, FILEA 

622 FORMAT (A10) 

OPEN(UNIT=21 ,ACCESS='APPEND' .FILE-FILEA) 

C 

C ONLY WRITE OUT MULTIPLES OF 5 DATA PAIRS 

C 

ITES=IMAX(2)/5 


n o 
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C 


999 

75 


IT-5*ITES 

TYPE *,IMAX(2) ,ITES,IT 
WRITE(21 ,*)IT,K2,K3,K4,P3 
APPEND TRANSMITENCES TO DATA FILE 

WRITE(21 ,999)( (ETCUR(I) .GAS (I) , J , J+4 ) ,J-1 ,IT,5) 
FORMAT(5(1X.F,1X,F)) 

CONTINUE 

RETURN 

END 


nonrionnonno nnnononnn no 
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APPENDIX B 
PROGRAM VFIT 


Program VFIT is listed in this appendix as discussed in 
the section Data Analysis. It is written to run in an 
interactive mode and uses plotting routines from the 
Tektronix PLOT10 graphics package and from Bevington [4]. 
The program is coded in Fortran- 10 and is designed to run on 
a 36 bit DecSystem-10 time shared computer under the TopslO 
monitor Version 6.03A. Some of the Format statements use 
the Fortran-10 free field characteristics ( F instead of 
Fw.d) and a few DO loops have indicies that are non-integer. 
Otherwise, this code conforms to the ANSI standard. 


PROGRAM VFIT 

DIMENSION KV(500) ,V(500) ,KABS(500) 
DIMENSION X(500) ,Y(500) ,SIGMAY(500) ,A(5) 
DIMENSION YFIT(300) 

DIMENSION ZVAL(300) ,ZIND(300) 

DIMENSION RES(20) ,YGUESS(20) 

REAL KV , KABS , LN2 , MASS , KBOLTZ , KO , KMAX 

REAL K1 ,K2 ,K3 

DOUBLE PRECISION XFILE 


THIS PROGRAM IS DESIGNED TO READ IN ABSORPTION 
COEFFICIENT DATA FOR A SINGLE LINE AND FIT 
A VOIGT FUNCTION TO IT. 


VOIGT FUNCION IS GENERATED BY THE FUNCTIONS 
K,K1,X2,K3 AS DESCRIBED BY B.H. ARMSTRONG 
IN JQSRT 7,PP 61-88. (1967) . 

THESE ARE ASSEMBLED IN THL FILE ARM. FOR 


WRITTEN BY DR. G.E. COPELAND AND L.N. MAJORANA 

DEPT OF PHYSICS 

OLD DOMINION UNIVERSITY 

NORFOLK, VA 23508 


DATA ARE IN FILE COF7.DAT 
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C INPUT DATA LOOKS LIKE THIS 

C 360, 12, 3, 79, 0.2600000 

C NPTS, MON, NDAY.NYEAR, PRESSURE IN TORR 

C 1129.3619 -0.0000081 1129.3623 -0.000016/ 1129.3626 

-0.0000303 1129.3629 -0.0000414 1129.3632 -0.0000426 

C FREQUENCY IN CM-1 FOLLOWED BY ABS COEFFICIENT BASE 

E 

C 5 DATA PAIRS PER LINE 

C 

C OPEN THE DATA FILE 

TYPE 82 

82 FORMATC INPUT THE NAME OF THE DATA 
FILE<COF? .DAT>? ’$) 

ACCEPT 83,XFILE 

83 FORMAT (A10) 

OPEN (UNIT-21 ,ACCESS-’SEQIN’ ,FILE-XFILE) 

READ(21 ,^)NPTS ,MON ,NDAY ,NYEAR,PRES 
TYPE 381 ,PRES,XFILE 

J81 FORMATC PRESSURE -’,F5.2,' TORR DATA FROM ’ ,A10) 

DO 10 1-1 , NPTS, 5 

READ(21 .*) (V(J) ,KABS(J) ,J-1 ,1+4) 

10 CONTINUE 

CLOSE(UNIT-21 ) 

TYPE 71 

71 FORMATC INPUT THE MODE FOR THE FILE?<1 ,0 , - 1 >? ’ $) 

ACCEPT *,MODE 
C 
C 

C FIND THE CENTER OF THE SPECTRAL LINE 

C 

YMAX— 1 0 . 

DO 20 1-1 ,NPTS 
IF(KABS(I) ,GT. YMAX)JJ=I 
IF(KABSd) .GT. YMAX)YMAX-KABS(I) 

20 CONTINUE 

TYPE 221 JJ YMAX 

221 FORMATC ’ MAX OCCURS AT ',13,' AND VALUE IS 

• ,1PE) 

C 

C 

C REMOVE 1129 CM-1 

DO 22 1=1 ,NPTS 
22 V(I)=V(I)-1 129. 

C 

C YMAX IS THE MAX OF THE ABS COEFFICIENT 

C JJ IS POINT IN DATA IT OCCURS 

C 

C TYPE *,JJ,YMAX 

C GRAB ABOUT 10 POINTS AROUND JJ AND FIT A QUADRATIC 

C TO THE DATA BY CALLING POLFIT . . . 

C Y=A(1 )+A(2'»*X+A(3)*X*X 
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C THEN THE POSITION OF THE MAX OF POLYNOMINAL 

C WILL BE THE CENTER LINE FOR ALL FUTURE CALCULATIONS 

C 

C GENERATE THE X AND Y ARRAYS FOR FITTING 

C 

NST-JJ-3 

NFIN-JJ+3 

NK-0 

DO 30 I-NST.NFIN 

NK-NK+1 

X(NK)-V(I) 

30 Y(NK)-KABSd) 

C 

C TYPE *,X,Y 

C PAUSE 

C 

C GENERATE THE MEAN AND ESTIMATED ERRORS FOR A SET OF 

DATA 

CALL XFIT(Y,SIGMAY,NK,MODE,YMEAN,SIGMAM,3IGM' 

DO 833 1-1 ,NK 

833 SIGMAY(I)-SQRT(( Y(I)-YMEAN)**2/(FLOA^ .^-1.)) 

C 

C BEVINGTON P76 

C 

C NOW CALL POLFIT BEVINGTON PI 40 

CALL POLFIT(X,Y,SIGMAY,NK,3 ,MODE,A,CHISQR) 

DO 49 1-1 ,NK 

ZVAL(I)-A(1 )+A(2)*X(I)+A(3)*X(I)**2 
C 1+A(4)*X(I)**3! +A(5)*X(I)**4 

49 CONTINUE 

C TYPE *,X,Y,ZVAL 

PAUSE 

CALL PLOT (X , Y , NK , X . ZVAL , 2 . PRFS , XFILE ) 
r PAii<^F 

C TYPE *,CHISQR,A 

C 

C CENTER OF THE LINE IS AT -A(2) / (2 .*A(3) ) 

VCEN— A(2)/(2.*A(3)) 

C FMAX-A(1 )+A(2)*VCEN+A(3)*VCEN**2 

TYPE 33.VCEN 

33 FORMATC CENTER FREQUENCY -’,F) 

C 

C STRIP OUT THE LINE CENTERS 

DO 50 1-1 ,NPTS 

50 V(I)-V(I)-VCEN 
C 

CALCULATE THE DOPPLER WIDTH 
PI-3.141592654 
SQPI-SQRT(PI) 

C 

r^.SS-3.*1 5 .9994 


C-2.99/9250E10 

C SPEED OF LIGHT IN CM/S 

KBOLTZ - 1.380622E-16 
C BOLTZMAN CONSTANT ERG/K 

T-285. 

C TEMP IN KELVIN 

AMU-1 .660531 E-24 
C 1AMU IN GRAMS 

LN2-ALOG(2.) 

SRLN2-SQRT(LN2) 

C 

AD-SQRT( 2 . *KBOLTZ*T*LN2 / (MASS*AMU) ) /C 
AD-AD*(1 129.+VCEN) 

AD1-AD*1000. 

AD2-AD1*C/1 .E-t-9 
CTYPE *,AD2 

TYPE 331 ADI AD2 

331 FORMAT (' DOPPLER WIDTH<FWHM> -'1PE* 

MILLIKAYSERS ' 0PF8 . 3 ’ MHZ ’ ) 

C 

CALCULATE THE 1ST APPROXIMATION TO THE LINE STRENGTH S 
C 

C NEED TO GET A DIFFERENT PIECE THAN WE DID FOR THE 

LINE CENTER 

C DETERMINATION THE NARROWER THE LINE THE LESS 

WE NEED 

C THE WIDER THE LINE THE MORE WE NEED. 

C THIS IS DETERMINED BY PRESSURE 

C WE COULD LOOK AT THE PRESSURE BUT INSTEAD LETS 

C COUNT THE NUMBER OF POINTS ABOVE 

C THE 10% LEVEL OF THE PEAK 

C IF THIS IS GREATER THAN 12 WE WILL GET 

C ONLY +/- 25 POINTS 

C 

C REPORT THE CUTOFF CRITERIAN CUTOFF VALUE 

C 

C 

C IF THIS IS LARGER WE WILL GET MORE POINTS FOR A 

BETTER 

C WING FIT 

CUTOFF-0 . 1 *YMAX 
TYPE 44 .CUTOFF 

44 FORMATC CUTOFF IS 10 PERCENT OF YMAX ’.F) 

ICUT-C 

CO 469 1-1 .NPTS 

IF( KABS(I) .GT. CUTOFF) ICUT-ICUT+1 
469 CONTINUE 

TYPE *,ICUT 
PAUSE 
I EXTRA- 30 


* ■i-n* 
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IF( ICUT.GE. 16) lEXTRA-AO 

NST-JJ-IEXTRA 

NFIN-JJ+IEXTRA 

NK-0 

DO 303 I-NST,NFIN 

NK-NK+1 

X(NK)-V(I) 

303 Y(NK)-KABSd) 

C 

C 

S-AREA(X,Y,NK,3) 

C TYPE *,S 

C 

C GENERATE KO 

C 

K0-(S/AD)*SRLN2/SQPI 
TYPE 990, KO 

990 FORMATC DIMENSIONAL CONSTANT OF THE REDUCED 

ABSORPTION 

1 COEFFICIENT , 1 PE) 

C 

C GENEPJVTE THE VARIBLE X 

C 

C CUT OFF POINTS ON FAR WINGS 

JL-JJ-IEXTRA 
JR-JJ+IEXTRA 
NP-0 

DO 88 I-JL.JR 
NP-NP+1 

88 X(NP)-V(I)*SRLN2/AD 

C TYPE *,V 

C 

TYPE 451 ,NP 

451 FORMATC TOTAL NUMBER OF DATA PAIRS USED IN 

ANALYSIS -'.U) 

C 

C GENERATE THE KABS DATA IN NORMAL FORM KABS/KO 

NP-0 

DO 89 I-JL,JR 
NP-NP+1 

89 KV (NP)-KABS(I)/KO 
C 

CALL A PLOTTING ROUTINE 

C CALL PLOT (X,KV,NP,X,KV,0, PRES, XFILE) 

C 

C 

C GUESS WHAT THE INTIAL VALUE OF Y =AL*SQLN2/AN IS??? 

r 

C USE THE FUNCTIONS K(X,Y) WITH X=0 .0 TO GUESS THE 

START OF Y 
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C 

C KMAX- VALUE OF THE K FUN 

KMAX— 1 0 . 

DO 41 1-1 ,NP 

IF(KMAX .LT. KV(I) ) KMAX - KV(I) 

41 CONTINUE 

TYPE 991 ,KMAX 

991 FORMAT(- MAX VALUE OF KABS/KO(AT LINE CENTER) 

-’.t) 

NZ-0 

DO 66 Z-0. 0,3. 0,0. 025 

NZ-NZ+1 

ZIND(NZ)-Z 

66 ZVAL(NZ)-ARM(0.0,Z) 

C TYPE *,ZVAL 

C 

CALL A ROUTINE TO INTERPOLATE TO GET CLOSER TO THE REAL Y 

VALUE 

C 

CALL INTER? ( ZVAL , ZIND , NZ , 3 . KMAX , YOUT) 

TYPE 99 2, YOUT 

992 FORMATC INTERPOLATED VALUE OF Y -'.F) 

DO 661 1-1 ,NP 

661 ZVAL(I)-ARM(ABS(X(D) ,YOUT) 

C TYPE *,ZVAL 

PAUSE 

CALL PLOT (X,KV,NP,X, ZVAL, 2 , PRES, XFILE) 

C 

C NOW WE HAVE A GOOD GUESS FOR THE VOIGT PROFILE 

C 

C WILL CHANGE THE Y VALUE OVER 10 VALUES OR SO AND 

CALC' 'LATE 

C THE CHI SQUARE TO THE DATA AND FIND MINIMUM Y FOR 

THOSE AND REDO 

C UNTIL CHI SQUARE IS MINIMUM 

C 

C X ARRARY HAS THE INDEPENDENT VARIABLES IN REDUCED 

FORM. 

C KV ARRAY IS THE REDUCED DATA 

C ZVAL IS THE FITTED DATA 

C NP IS NUMBER OF DATA POINTS 

C 

C 1ST START UP THE CHI SQR CALCULATION USING YOUT AS 

THE 

C 1ST VALUE OF THE K(X , Y) ---THEORY FIT ABOVE 

DO 551 1-1 ,NP 
551 SIGMAY(I)-0.0 

C^LL XFIT(KV , SIGMAY ,NP .MODE , YMEAN , SIGMAM, SIGMA) 

C 

CALCULATE THE SIGMAY'S 
C 
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DO 561 1-1 ,NP 

561 SIGMAY(I)-SQRT( ( KV(I)-YMEAN)**2/(FL0AT(NP)-1 .) 

) 

NFREE-NP-2 

REST-FCHISQ (KV , SIGMAY , NP . NFREE , MODE , ZVAL) 

C 

C REST IS THE REDUCED CHI SQUARE FOR THE GUESS 

:ST-YOUT-YOUT/10. 

YEND*«YOUT-*-YOUT/10. 

YSTEP-YOUT/100. 

95 I-O 

DO 552 YS-YST,YEND,YSTEP 
I-I+1 

552 YGUESS(I)-YS 
ITOTAL -I 

I-O 

C 

CALL THE FCHISQ ROUTINE 

TYPE 441 ,I,REST,YOUT 
DO 553 1-1, ITOTAL 

C 

CALCULATE THE ZVAL'S FOR THIS YGUESS 
DO 554 J-1 ,NP 

554 ZVAL(J)-ARM( ABS(X(J)) .YGUESS(I) ) 

C 

C GENERATE THE SIGMAY 'S 

DO 5551 IK-1 ,NP 
5551 SIGMAY(IK)-0 .0 

DO 5661 IH-1 ,NP 

5661 SIGMAY(IH)-SQRT( ( 

KV(IH)-YMEAN)**2/(FLOAT(NP)-1 .) ) 

CALL XFIT(KV , SIGMAY ,NP .MODE ,YMEAN , SIGMAM, SIGMA) 

RES (I) -FCHISQ (KV , SIGMAY , NP , NFREE , MODE , ZVAL) 

TYPE 441 ,I ,RES(I) , YGUESS (I) 

441 FORMATC TRY ’,12,' REDUCED CHISR-'IPE’Y ='1PE) 

C 

C 

C SKIP OUT IF REDUCED CHISQ DOES NOT CHANGE BY 0 .C05 

C IF( I.EQ.1) GOTO 553 

C 
C 

C IF(ABS( ( RES(I)-RES(I-1 ) )/RES(I) ) .LT. 0.005) 

IKEEP-I 

C IF(ABS( ( RES(I)-RES(I-1 ) )/RES(I) ) .LT 0.01) 

GOTO 6812 

553 CONTINUE 

68 CONTINUE 

YMIN-10. 

DO 69 1=1 ITOTAL 

IF( RES(I) .LT. YMIN) IKEEP-I 

IF( RES(I) .LT. YMIN) YMIN=RES(Iy 


69 CONTINUE 

IF( IKEEP.EQ. ITOTAL) T\PE 510 
510 FORMATC OOPS! STILL DECREASING AT EDGE OF 

SEARCH. I WILL 

1 CONTINUE ') 

IF( IKEEP.EQ. ITOTAL) VST-YEND 

IF( IKEEP.EQ. ITOTAL) YEND-YST+YOUT/ 1 0 . 

IF( IKEEP.EQ. ITOTADGOTO 95 
6812 TYPE 690,YGUESS(IKEEP) 

690 FORMATC RATIO OF LORENTZ TO DOPPLER WIDTH 
-• ,1PE) 

AL- YGUESS ( IKEEP ) *AD / S RLN2 
AL1-AL*1000. 

AL2-AL1*C/1 .E-»-9 
TYPE 691 ,AL1 ,AL2 

691 FORMATC LORENTZ WIDTH '.IPE,' MILLIKAYSERS 
-•0PF8.3,' MHZ’) 

C PAUSE 

C 

CPREPARE A RICHER ARRAY FOR PLOTTING MAKE IT 3*NP LONG 
C 

C 1ST DO THE X \RRAY 

C PUT NEW X ARRAY IN V 

C PUT DATA Y ARRAY IN KABS 

C PUT NEW FITTED DATA IN Y 

MP3-3*NP 
DO 46 1-1 .NPTS 

46 KABS(I)-0.0 
DO 461 1-1 ,NPTS 

461 V(I)-0.0 

IL-1 

DO 4/ J-1 ,NP 
V(IL)-X(J) 

V(IL+1 )-X(J) 

V(IL+2)-X(J) 

KABSdL)-KV(J) 

KABS(IL+1 )-KV(J) 

KABS(IL+2)-Kv(J) 

IL-IL+3 

47 CONTINUE 

C TYPE *,V,KABS 

C 

CALCULATE A RIC^ NEW X ARRY FOR INDEPENDENT VAR^.ABLE 
C FOR THE SECOND GRA.PH OF FITTED VS. X 

C 

XMIN-V( I ) 

XMAX=V(NP3) 

C 

C STOE NEW INDEPENT VARIABLE IN Y 

C 
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XINC-2 . ★XMAX/ FLOAT (NP3 ) 

IC-0 

DO 201 XO-XMIN,XMAX,XINC 
IC-IC+I 

20 1 Y(IC)-XO 

C TYPE *,IC,NP3 

C PAUSE 

C 

CALL ROUTINE TO GET THE NEW DEPENDENT VARIBLE AT THE SPACING 
C DESIRED XINC 

n 

C STORE IN ZVAL 

DO 345 1-1 ,NP3 

345 ZIND(I)-ARM( ABS(Y(I)) , YGUESS(IKEEP)) 

C 

PAUSE 

CALL PLOT(V,KABS,NP3,y,ZIND,2,PRES,XFILE) 

END 

C 

C - 

C 

SUBROUTINE PLOT(X . Y ,NP ,X1 ,Y1 , IPLT , P .XFILE) 

DIMENSION X(1),Y(1) 

DIMENSION XI (1 ) ,Y1 (1 ) 

DIMENSION YLABEL(18) ,IYLAB(18) 

DIMENSION TR(4) 

DIMENSION XLABEL(48) ,IXLAB(48) 

DIMENSION TITLE(22) ,ITIT(22) 

DIMENSION TITLE2(28) ,ITIT2(28) 

DOUBLE PRECISION XFILE 
C 


C 

C 


• , 'F' 
C22 


DATA FOR THE TITLES FOR THE GRAPHS 


DATA TITLE/ ’V* , ’0 
, 'U’ , 'N' , 'C , 'T' , 'I • 
1 , ' • , 'F' , 'I ■ , 'T* 
CHARACTERS 



'/ 


DATA TR/'T' , '0 ‘ ’ R ' , ' R ‘ / 

C 4 CHARACTERS 

DATA 

XLABEL/ ’W ,'A’,'V',’EV, 'N’,’U‘,'M','B’,'E’,’R' 


! 'S' . 'C , 'A' , 'L’ , 'E' , ' [ ■ , ’D’ , 'O- , ’P' 


■L- ’E’ 'R‘ 


2'W','I','D’,'T','H','=','0','.'.'9',’8','5','M','K'.’=',*2’, 
3'9’ , • . ' , '5' , 'M' , 'H- , 'Z' . ■ 1 •/ 

C48 CHARACTERS 

DATA YLABEL/ 'R-,'E','D','U','C','E’,’D’,' 

, A , 5 , b , . , 

1 'C , 'O’ , 'E' , 'F' , 'F' , ' . ' / 

Cl 8 CHARACTERS 
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DATA TITLE2/'0* , '3' 

• , 'A' , ’T* , '1 • , '1 • , 

I '2* , •9' , • . ' '4' , '4 

• , *F\ 'R', 'O’ , 'M' , ' V 
C28 CHARACTERS 

CALL IN ITT (240) 

CALL BINITT 


' , 'L’ , , ’N* , 'E' , • 


3' , *C’ , 'M* , , ’T , • 


C 

C 


CALL NPTS(NP) 

CALL XFRM(2) 

CALL Y?RM(2) 

IFCIPLT .NE.O) CALL 
IF(I’'LT .NE.O) CALL 
CALL CHECK(X.Y) 

CALL DSPLAY(X.Y) 
IF(IPLT .EQ.2) CALL 
IF(IPLT .EQ.2) CALL 
IF( IPLT.EQ.2) CALL 


SYMBLd ) 
LINE(-I) 


SYIIBLCO) 
LINE(O) 
CPLOT(X1 ,Y1 ) 


MOVE AND PUT ON LABELS 


CALL MOVABS(20 ,600) 

CALL KA1 2AS(18 .YLABEL.IYLAB) 
CALL VLABEL(18 .lYLAB) 

CALL KA12AS(48,XLABEL,IXLAB) 
CALL NOTATE(325 ,20,48, IXLAB) 
775 CONTINUE 

CALL KA12AS(28,TITLE2,ITIT2) 
CALL MOVABS(200,750) 

CALL HLABEL(28,ITIT2) 

CALL ANMODE 
TYPE 10,XFILE 

10 FORMAT(1X,A10) 

CALL KA12AS (22, TITLE, ITIT) 
CALL MOVABS(200 ,725) 

CALL HLABEL(22,ITIT) 

CALL ANMODE 
TYPE 1 1 ,P 

11 FORMAT(1X,F8.2) 

CALL BELL 

CALL TINPUT(IEND) 

CALL ERASE 


RETURN 

END 


C-- 


CODE TO CALCULATE THE VOICT PROFILE 
C 


C TAKEN FROM 

C SPECTRUM LINE PROFIESrTHE VOIGT FUNCTION 

C BY B.H. ARMSTRONG 

C JQSRT 7 ,PP61-88 , (1967) 

C 
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FUNCTION ARM(X.Y) 

C THIS FUNCTION IS THF REAL PART OF THE COMPLEX 

PROBILITY FUNCTION 

C OR THE VOIGT SPECTRUM LINE PROFILE 

REAL K1 ,K2,K3 
COMMON W(10) ,T(10) ,Y2 

DATA W/4.62243670E-1 .2.86675505E-1 ,1 .^^0 I 7206E-1 , 

1 2.48105209E-2,3.24377334E-3,2.283386j6E-4 , 

2 

7 .80255648E-6.1 .08606937E-7 ,4 .39934099E- 10 ,2 .22939365E-1 3/ 
DATA 

T/0 .245340708,0.737473729,1 .23407622,1 .73853771 ,2.25497400, 
1 

2 . /:i880606 ,3 .347P5457 ,3 .94476404 ,4 .60368245 ,5 .38748089/ 
Y2-Y*Y 
C 

C TESTS 

C 

IF(Y.LT.1 .0.AND.X.LT.4 .0 .OR.Y.LT.1 .8/(X+1 .0)) GOTO 

300 

IF(Y.LT.2 .5.AND.X.LT.4 .0) GOTO 200 
100 ARM-K3(X,Y) 

RETURN 

200 ARM-K2(X,Y) 

RETURN 

300 ARM-K1(X,Y) 

RETURN 

END 

C 

C 

C 

FUNCTION K1(X,Y) 

REAL K1 

DOUBLE PRECISION C (34 ) , COEF , BNO 1 , BN02 , BN , XI , F 
F3 (T) -EXP (T**2-X**2 ) 

C FROM HERE TO STATEMENT 30 WE CALCULATE DAWSONS 

FUNCTION 

C ENTER HUMMERS CHEBYSHEV COEFFICIENTS C(I) 

DATA C/ .1999999999972224 ,-.1840000000029998 , 

1 

. 1 558399999965025 ,- . 1216640000043988 , .0877081599940391 , 

2 

-.0585141248086907, .03621 5730 16239 14, -.020849 7654398036 , 

3 

.01 1 19601 16346270, 56231 896 1671 09D-2, .26487634 172265D-2, 

4 

-.1 1 732670757 704D-2, .48995 199 78088D-3,-. 193363080 1528D-3 , 

5 

. 72287 /44678PD-4 ,-.256555 124979D-4 , .86620 73684 1 D-5 , 

6 - . 278763797 1 9D-5, .8566873627D-6,-. 25 18433784D-6, 
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7 

. 709360221 D-/ ,-.1917322570-7, .49801 256D-8,-. 1244/734 D-8, 

8 

. 299 / 7 / 7D-9 , - .6964500- 1 0 ,. 1 562620- 1 0 , - . 3389 /O- 1 1 ,. 7 1 1 60- 1 2 , 

9 -.14470-12, .2850- 13, -.550- 14, .100- 14, -.20-1 5/ 
Y2-Y*Y 

IF((X**2-Y2) .GT.70.0)GO TO 2 
U1 -EXP(-X**2+Y2)*COS(2 .*X*Y) 

GO TO 5 
2 Ul-0.0 

5 IF(X.GT.5.0) GO TO 1000 

C CLENSKAWS ALGORITHM AS GIVEN BY HUMMER 

BN01-0.000 
BN02-0.000 
Xl-X/5.000 

COEF-4 .000*X1 **2-2.000 

00 20 1-1 ,34 
II-35-I 

BN-COEF*BN01-BN02-*-C(II) 

BN02-BN01 
20 BN01-BN 

30 F-X1*(BN-BN02) 

40 ONOl-1 .0-2.0*X*SNGL(F) 

1100 0N02-SNGL(F) 

GO TO 1200 

1000 

0N01— (.5/X**2+.75/X**4+1 .875/X**6+6 .5625/X**8+29 .531 25/X**1 0 

•f 

1 162. 4218/X**12+1055. 7421 /X**14) 

0N02-(1 .-ON01 )/(2.*X) 

1200 FUNCT-Y*DN01 

IF(Y.LE.1 .OE-8) GO TO 2500 

Q-1 .0 

YN-Y 

DO 2000 1-2,50 

DN-(X*DN01+DN02)*(-2.)/FLOAT(I) 

DN02-DN01 

DN01-DN 

IF (MOD(I,2))2000, 2000, 1500 
1 500 Q— Q 

YN-YN*Y2 

G-DN*YN 

FUNCT-FUNCT+Q*G 

IF (ABS(G/FUNCT) .LE.l .00E-08)GO TO 2500 
2000 CONTINUE 

2500 K1-U1-1 .1283791 7*FUNCT 

RETURN 
END 
C 

C - 

C 


o on 
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FUNCTION K2(X,Y) 

RFAL K2 

COMMON W(10) ,T(10) ,Y2 
G-0.0 

DO 100 1-1,10 


100 

G-G+ (4 . *T ( I ) **2 - 2 . ) * ( R*ATAN ( R / Y ) +S*ATAN ( S / Y ) - . 5 *Y* ( ALOG ( Y 2 +R* 
★ 2 ) 

1 +ALOG(Y2+S**2)))*W(I) 

K2-0.318309886*G 

RETURN 

END 


FUNCTION K3(X,Y) 

REAL K3 

COMMON W(10) ,T(10) ,Y2 
G-0.0 

DO 100 1=1,10 

100 

G-G+(1 .0E0/((X-T(I))**2+Y2)+1 .OEO / ( (X+T(I ) )**2+Y2 ) )*W(I ) 
K3-0.318309886*Y’^G 
RETURN 
END 




